|
kusano |
7d535a |
/*! @file get_perm_c.c
|
|
kusano |
7d535a |
* \brief Matrix permutation operations
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* -- SuperLU routine (version 3.1) --
|
|
kusano |
7d535a |
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
|
|
kusano |
7d535a |
* and Lawrence Berkeley National Lab.
|
|
kusano |
7d535a |
* August 1, 2008
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
#include "slu_ddefs.h"
|
|
kusano |
7d535a |
#include "colamd.h"
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
extern int genmmd_(int *, int *, int *, int *, int *, int *, int *,
|
|
kusano |
7d535a |
int *, int *, int *, int *, int *);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
void
|
|
kusano |
7d535a |
get_colamd(
|
|
kusano |
7d535a |
const int m, /* number of rows in matrix A. */
|
|
kusano |
7d535a |
const int n, /* number of columns in matrix A. */
|
|
kusano |
7d535a |
const int nnz,/* number of nonzeros in matrix A. */
|
|
kusano |
7d535a |
int *colptr, /* column pointer of size n+1 for matrix A. */
|
|
kusano |
7d535a |
int *rowind, /* row indices of size nz for matrix A. */
|
|
kusano |
7d535a |
int *perm_c /* out - the column permutation vector. */
|
|
kusano |
7d535a |
)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
int Alen, *A, i, info, *p;
|
|
kusano |
7d535a |
double knobs[COLAMD_KNOBS];
|
|
kusano |
7d535a |
int stats[COLAMD_STATS];
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Alen = colamd_recommended(nnz, m, n);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
colamd_set_defaults(knobs);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if (!(A = (int *) SUPERLU_MALLOC(Alen * sizeof(int))) )
|
|
kusano |
7d535a |
ABORT("Malloc fails for A[]");
|
|
kusano |
7d535a |
if (!(p = (int *) SUPERLU_MALLOC((n+1) * sizeof(int))) )
|
|
kusano |
7d535a |
ABORT("Malloc fails for p[]");
|
|
kusano |
7d535a |
for (i = 0; i <= n; ++i) p[i] = colptr[i];
|
|
kusano |
7d535a |
for (i = 0; i < nnz; ++i) A[i] = rowind[i];
|
|
kusano |
7d535a |
info = colamd(m, n, Alen, A, p, knobs, stats);
|
|
kusano |
7d535a |
if ( info == FALSE ) ABORT("COLAMD failed");
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) perm_c[p[i]] = i;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SUPERLU_FREE(A);
|
|
kusano |
7d535a |
SUPERLU_FREE(p);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
/*! \brief
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Purpose
|
|
kusano |
7d535a |
* =======
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Form the structure of A'*A. A is an m-by-n matrix in column oriented
|
|
kusano |
7d535a |
* format represented by (colptr, rowind). The output A'*A is in column
|
|
kusano |
7d535a |
* oriented format (symmetrically, also row oriented), represented by
|
|
kusano |
7d535a |
* (ata_colptr, ata_rowind).
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* This routine is modified from GETATA routine by Tim Davis.
|
|
kusano |
7d535a |
* The complexity of this algorithm is: SUM_{i=1,m} r(i)^2,
|
|
kusano |
7d535a |
* i.e., the sum of the square of the row counts.
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Questions
|
|
kusano |
7d535a |
* =========
|
|
kusano |
7d535a |
* o Do I need to withhold the *dense* rows?
|
|
kusano |
7d535a |
* o How do I know the number of nonzeros in A'*A?
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
void
|
|
kusano |
7d535a |
getata(
|
|
kusano |
7d535a |
const int m, /* number of rows in matrix A. */
|
|
kusano |
7d535a |
const int n, /* number of columns in matrix A. */
|
|
kusano |
7d535a |
const int nz, /* number of nonzeros in matrix A */
|
|
kusano |
7d535a |
int *colptr, /* column pointer of size n+1 for matrix A. */
|
|
kusano |
7d535a |
int *rowind, /* row indices of size nz for matrix A. */
|
|
kusano |
7d535a |
int *atanz, /* out - on exit, returns the actual number of
|
|
kusano |
7d535a |
nonzeros in matrix A'*A. */
|
|
kusano |
7d535a |
int **ata_colptr, /* out - size n+1 */
|
|
kusano |
7d535a |
int **ata_rowind /* out - size *atanz */
|
|
kusano |
7d535a |
)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
register int i, j, k, col, num_nz, ti, trow;
|
|
kusano |
7d535a |
int *marker, *b_colptr, *b_rowind;
|
|
kusano |
7d535a |
int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if ( !(marker = (int*) SUPERLU_MALLOC((SUPERLU_MAX(m,n)+1)*sizeof(int))) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for marker[]");
|
|
kusano |
7d535a |
if ( !(t_colptr = (int*) SUPERLU_MALLOC((m+1) * sizeof(int))) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC t_colptr[]");
|
|
kusano |
7d535a |
if ( !(t_rowind = (int*) SUPERLU_MALLOC(nz * sizeof(int))) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for t_rowind[]");
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Get counts of each column of T, and set up column pointers */
|
|
kusano |
7d535a |
for (i = 0; i < m; ++i) marker[i] = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i)
|
|
kusano |
7d535a |
++marker[rowind[i]];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
t_colptr[0] = 0;
|
|
kusano |
7d535a |
for (i = 0; i < m; ++i) {
|
|
kusano |
7d535a |
t_colptr[i+1] = t_colptr[i] + marker[i];
|
|
kusano |
7d535a |
marker[i] = t_colptr[i];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Transpose the matrix from A to T */
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j)
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
col = rowind[i];
|
|
kusano |
7d535a |
t_rowind[marker[col]] = j;
|
|
kusano |
7d535a |
++marker[col];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* ----------------------------------------------------------------
|
|
kusano |
7d535a |
compute B = T * A, where column j of B is:
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Struct (B_*j) = UNION ( Struct (T_*k) )
|
|
kusano |
7d535a |
A_kj != 0
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
do not include the diagonal entry
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
( Partition A as: A = (A_*1, ..., A_*n)
|
|
kusano |
7d535a |
Then B = T * A = (T * A_*1, ..., T * A_*n), where
|
|
kusano |
7d535a |
T * A_*j = (T_*1, ..., T_*m) * A_*j. )
|
|
kusano |
7d535a |
---------------------------------------------------------------- */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Zero the diagonal flag */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) marker[i] = -1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* First pass determines number of nonzeros in B */
|
|
kusano |
7d535a |
num_nz = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
/* Flag the diagonal so it's not included in the B matrix */
|
|
kusano |
7d535a |
marker[j] = j;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
/* A_kj is nonzero, add pattern of column T_*k to B_*j */
|
|
kusano |
7d535a |
k = rowind[i];
|
|
kusano |
7d535a |
for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
|
|
kusano |
7d535a |
trow = t_rowind[ti];
|
|
kusano |
7d535a |
if ( marker[trow] != j ) {
|
|
kusano |
7d535a |
marker[trow] = j;
|
|
kusano |
7d535a |
num_nz++;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
*atanz = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Allocate storage for A'*A */
|
|
kusano |
7d535a |
if ( !(*ata_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for ata_colptr[]");
|
|
kusano |
7d535a |
if ( *atanz ) {
|
|
kusano |
7d535a |
if ( !(*ata_rowind = (int*) SUPERLU_MALLOC( *atanz * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for ata_rowind[]");
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
b_colptr = *ata_colptr; /* aliasing */
|
|
kusano |
7d535a |
b_rowind = *ata_rowind;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Zero the diagonal flag */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) marker[i] = -1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Compute each column of B, one at a time */
|
|
kusano |
7d535a |
num_nz = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
b_colptr[j] = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Flag the diagonal so it's not included in the B matrix */
|
|
kusano |
7d535a |
marker[j] = j;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
/* A_kj is nonzero, add pattern of column T_*k to B_*j */
|
|
kusano |
7d535a |
k = rowind[i];
|
|
kusano |
7d535a |
for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
|
|
kusano |
7d535a |
trow = t_rowind[ti];
|
|
kusano |
7d535a |
if ( marker[trow] != j ) {
|
|
kusano |
7d535a |
marker[trow] = j;
|
|
kusano |
7d535a |
b_rowind[num_nz++] = trow;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
b_colptr[n] = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SUPERLU_FREE(marker);
|
|
kusano |
7d535a |
SUPERLU_FREE(t_colptr);
|
|
kusano |
7d535a |
SUPERLU_FREE(t_rowind);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/*! \brief
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Purpose
|
|
kusano |
7d535a |
* =======
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Form the structure of A'+A. A is an n-by-n matrix in column oriented
|
|
kusano |
7d535a |
* format represented by (colptr, rowind). The output A'+A is in column
|
|
kusano |
7d535a |
* oriented format (symmetrically, also row oriented), represented by
|
|
kusano |
7d535a |
* (b_colptr, b_rowind).
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
void
|
|
kusano |
7d535a |
at_plus_a(
|
|
kusano |
7d535a |
const int n, /* number of columns in matrix A. */
|
|
kusano |
7d535a |
const int nz, /* number of nonzeros in matrix A */
|
|
kusano |
7d535a |
int *colptr, /* column pointer of size n+1 for matrix A. */
|
|
kusano |
7d535a |
int *rowind, /* row indices of size nz for matrix A. */
|
|
kusano |
7d535a |
int *bnz, /* out - on exit, returns the actual number of
|
|
kusano |
7d535a |
nonzeros in matrix A'*A. */
|
|
kusano |
7d535a |
int **b_colptr, /* out - size n+1 */
|
|
kusano |
7d535a |
int **b_rowind /* out - size *bnz */
|
|
kusano |
7d535a |
)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
register int i, j, k, col, num_nz;
|
|
kusano |
7d535a |
int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
|
|
kusano |
7d535a |
int *marker;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if ( !(marker = (int*) SUPERLU_MALLOC( n * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for marker[]");
|
|
kusano |
7d535a |
if ( !(t_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for t_colptr[]");
|
|
kusano |
7d535a |
if ( !(t_rowind = (int*) SUPERLU_MALLOC( nz * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails t_rowind[]");
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Get counts of each column of T, and set up column pointers */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) marker[i] = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i)
|
|
kusano |
7d535a |
++marker[rowind[i]];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
t_colptr[0] = 0;
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) {
|
|
kusano |
7d535a |
t_colptr[i+1] = t_colptr[i] + marker[i];
|
|
kusano |
7d535a |
marker[i] = t_colptr[i];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Transpose the matrix from A to T */
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j)
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
col = rowind[i];
|
|
kusano |
7d535a |
t_rowind[marker[col]] = j;
|
|
kusano |
7d535a |
++marker[col];
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* ----------------------------------------------------------------
|
|
kusano |
7d535a |
compute B = A + T, where column j of B is:
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
Struct (B_*j) = Struct (A_*k) UNION Struct (T_*k)
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
do not include the diagonal entry
|
|
kusano |
7d535a |
---------------------------------------------------------------- */
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Zero the diagonal flag */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) marker[i] = -1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* First pass determines number of nonzeros in B */
|
|
kusano |
7d535a |
num_nz = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
/* Flag the diagonal so it's not included in the B matrix */
|
|
kusano |
7d535a |
marker[j] = j;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Add pattern of column A_*k to B_*j */
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
k = rowind[i];
|
|
kusano |
7d535a |
if ( marker[k] != j ) {
|
|
kusano |
7d535a |
marker[k] = j;
|
|
kusano |
7d535a |
++num_nz;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Add pattern of column T_*k to B_*j */
|
|
kusano |
7d535a |
for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
k = t_rowind[i];
|
|
kusano |
7d535a |
if ( marker[k] != j ) {
|
|
kusano |
7d535a |
marker[k] = j;
|
|
kusano |
7d535a |
++num_nz;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
*bnz = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Allocate storage for A+A' */
|
|
kusano |
7d535a |
if ( !(*b_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for b_colptr[]");
|
|
kusano |
7d535a |
if ( *bnz) {
|
|
kusano |
7d535a |
if ( !(*b_rowind = (int*) SUPERLU_MALLOC( *bnz * sizeof(int)) ) )
|
|
kusano |
7d535a |
ABORT("SUPERLU_MALLOC fails for b_rowind[]");
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Zero the diagonal flag */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) marker[i] = -1;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Compute each column of B, one at a time */
|
|
kusano |
7d535a |
num_nz = 0;
|
|
kusano |
7d535a |
for (j = 0; j < n; ++j) {
|
|
kusano |
7d535a |
(*b_colptr)[j] = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Flag the diagonal so it's not included in the B matrix */
|
|
kusano |
7d535a |
marker[j] = j;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Add pattern of column A_*k to B_*j */
|
|
kusano |
7d535a |
for (i = colptr[j]; i < colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
k = rowind[i];
|
|
kusano |
7d535a |
if ( marker[k] != j ) {
|
|
kusano |
7d535a |
marker[k] = j;
|
|
kusano |
7d535a |
(*b_rowind)[num_nz++] = k;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Add pattern of column T_*k to B_*j */
|
|
kusano |
7d535a |
for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
|
|
kusano |
7d535a |
k = t_rowind[i];
|
|
kusano |
7d535a |
if ( marker[k] != j ) {
|
|
kusano |
7d535a |
marker[k] = j;
|
|
kusano |
7d535a |
(*b_rowind)[num_nz++] = k;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
(*b_colptr)[n] = num_nz;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SUPERLU_FREE(marker);
|
|
kusano |
7d535a |
SUPERLU_FREE(t_colptr);
|
|
kusano |
7d535a |
SUPERLU_FREE(t_rowind);
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/*! \brief
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Purpose
|
|
kusano |
7d535a |
* =======
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* GET_PERM_C obtains a permutation matrix Pc, by applying the multiple
|
|
kusano |
7d535a |
* minimum degree ordering code by Joseph Liu to matrix A'*A or A+A'.
|
|
kusano |
7d535a |
* or using approximate minimum degree column ordering by Davis et. al.
|
|
kusano |
7d535a |
* The LU factorization of A*Pc tends to have less fill than the LU
|
|
kusano |
7d535a |
* factorization of A.
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* Arguments
|
|
kusano |
7d535a |
* =========
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* ispec (input) int
|
|
kusano |
7d535a |
* Specifies the type of column ordering to reduce fill:
|
|
kusano |
7d535a |
* = 1: minimum degree on the structure of A^T * A
|
|
kusano |
7d535a |
* = 2: minimum degree on the structure of A^T + A
|
|
kusano |
7d535a |
* = 3: approximate minimum degree for unsymmetric matrices
|
|
kusano |
7d535a |
* If ispec == 0, the natural ordering (i.e., Pc = I) is returned.
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* A (input) SuperMatrix*
|
|
kusano |
7d535a |
* Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
|
|
kusano |
7d535a |
* of the linear equations is A->nrow. Currently, the type of A
|
|
kusano |
7d535a |
* can be: Stype = NC; Dtype = _D; Mtype = GE. In the future,
|
|
kusano |
7d535a |
* more general A can be handled.
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
* perm_c (output) int*
|
|
kusano |
7d535a |
* Column permutation vector of size A->ncol, which defines the
|
|
kusano |
7d535a |
* permutation matrix Pc; perm_c[i] = j means column i of A is
|
|
kusano |
7d535a |
* in position j in A*Pc.
|
|
kusano |
7d535a |
*
|
|
kusano |
7d535a |
*/
|
|
kusano |
7d535a |
void
|
|
kusano |
7d535a |
get_perm_c(int ispec, SuperMatrix *A, int *perm_c)
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
NCformat *Astore = A->Store;
|
|
kusano |
7d535a |
int m, n, bnz = 0, *b_colptr, i;
|
|
kusano |
7d535a |
int delta, maxint, nofsub, *invp;
|
|
kusano |
7d535a |
int *b_rowind, *dhead, *qsize, *llist, *marker;
|
|
kusano |
7d535a |
double t, SuperLU_timer_();
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
m = A->nrow;
|
|
kusano |
7d535a |
n = A->ncol;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
t = SuperLU_timer_();
|
|
kusano |
7d535a |
switch ( ispec ) {
|
|
kusano |
7d535a |
case (NATURAL): /* Natural ordering */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) perm_c[i] = i;
|
|
kusano |
7d535a |
#if ( PRNTlevel>=1 )
|
|
kusano |
7d535a |
printf("Use natural column ordering.\n");
|
|
kusano |
7d535a |
#endif
|
|
kusano |
7d535a |
return;
|
|
kusano |
7d535a |
case (MMD_ATA): /* Minimum degree ordering on A'*A */
|
|
kusano |
7d535a |
getata(m, n, Astore->nnz, Astore->colptr, Astore->rowind,
|
|
kusano |
7d535a |
&bnz, &b_colptr, &b_rowind);
|
|
kusano |
7d535a |
#if ( PRNTlevel>=1 )
|
|
kusano |
7d535a |
printf("Use minimum degree ordering on A'*A.\n");
|
|
kusano |
7d535a |
#endif
|
|
kusano |
7d535a |
t = SuperLU_timer_() - t;
|
|
kusano |
7d535a |
/*printf("Form A'*A time = %8.3f\n", t);*/
|
|
kusano |
7d535a |
break;
|
|
kusano |
7d535a |
case (MMD_AT_PLUS_A): /* Minimum degree ordering on A'+A */
|
|
kusano |
7d535a |
if ( m != n ) ABORT("Matrix is not square");
|
|
kusano |
7d535a |
at_plus_a(n, Astore->nnz, Astore->colptr, Astore->rowind,
|
|
kusano |
7d535a |
&bnz, &b_colptr, &b_rowind);
|
|
kusano |
7d535a |
#if ( PRNTlevel>=1 )
|
|
kusano |
7d535a |
printf("Use minimum degree ordering on A'+A.\n");
|
|
kusano |
7d535a |
#endif
|
|
kusano |
7d535a |
t = SuperLU_timer_() - t;
|
|
kusano |
7d535a |
/*printf("Form A'+A time = %8.3f\n", t);*/
|
|
kusano |
7d535a |
break;
|
|
kusano |
7d535a |
case (COLAMD): /* Approximate minimum degree column ordering. */
|
|
kusano |
7d535a |
get_colamd(m, n, Astore->nnz, Astore->colptr, Astore->rowind, perm_c);
|
|
kusano |
7d535a |
#if ( PRNTlevel>=1 )
|
|
kusano |
7d535a |
printf(".. Use approximate minimum degree column ordering.\n");
|
|
kusano |
7d535a |
#endif
|
|
kusano |
7d535a |
return;
|
|
kusano |
7d535a |
default:
|
|
kusano |
7d535a |
ABORT("Invalid ISPEC");
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
if ( bnz != 0 ) {
|
|
kusano |
7d535a |
t = SuperLU_timer_();
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Initialize and allocate storage for GENMMD. */
|
|
kusano |
7d535a |
delta = 0; /* DELTA is a parameter to allow the choice of nodes
|
|
kusano |
7d535a |
whose degree <= min-degree + DELTA. */
|
|
kusano |
7d535a |
maxint = 2147483647; /* 2**31 - 1 */
|
|
kusano |
7d535a |
invp = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
|
|
kusano |
7d535a |
if ( !invp ) ABORT("SUPERLU_MALLOC fails for invp.");
|
|
kusano |
7d535a |
dhead = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
|
|
kusano |
7d535a |
if ( !dhead ) ABORT("SUPERLU_MALLOC fails for dhead.");
|
|
kusano |
7d535a |
qsize = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
|
|
kusano |
7d535a |
if ( !qsize ) ABORT("SUPERLU_MALLOC fails for qsize.");
|
|
kusano |
7d535a |
llist = (int *) SUPERLU_MALLOC(n*sizeof(int));
|
|
kusano |
7d535a |
if ( !llist ) ABORT("SUPERLU_MALLOC fails for llist.");
|
|
kusano |
7d535a |
marker = (int *) SUPERLU_MALLOC(n*sizeof(int));
|
|
kusano |
7d535a |
if ( !marker ) ABORT("SUPERLU_MALLOC fails for marker.");
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Transform adjacency list into 1-based indexing required by GENMMD.*/
|
|
kusano |
7d535a |
for (i = 0; i <= n; ++i) ++b_colptr[i];
|
|
kusano |
7d535a |
for (i = 0; i < bnz; ++i) ++b_rowind[i];
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
genmmd_(&n, b_colptr, b_rowind, perm_c, invp, &delta, dhead,
|
|
kusano |
7d535a |
qsize, llist, marker, &maxint, &nofsub);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Transform perm_c into 0-based indexing. */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) --perm_c[i];
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SUPERLU_FREE(invp);
|
|
kusano |
7d535a |
SUPERLU_FREE(dhead);
|
|
kusano |
7d535a |
SUPERLU_FREE(qsize);
|
|
kusano |
7d535a |
SUPERLU_FREE(llist);
|
|
kusano |
7d535a |
SUPERLU_FREE(marker);
|
|
kusano |
7d535a |
SUPERLU_FREE(b_rowind);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
t = SuperLU_timer_() - t;
|
|
kusano |
7d535a |
/* printf("call GENMMD time = %8.3f\n", t);*/
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
} else { /* Empty adjacency structure */
|
|
kusano |
7d535a |
for (i = 0; i < n; ++i) perm_c[i] = i;
|
|
kusano |
7d535a |
}
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
SUPERLU_FREE(b_colptr);
|
|
kusano |
7d535a |
}
|