kusano 7d535a
kusano 7d535a
/*! @file 
kusano 7d535a
 * \brief Finds a row permutation so that the matrix has large entries on the diagonal
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 4.0) --
kusano 7d535a
 * Lawrence Berkeley National Laboratory.
kusano 7d535a
 * June 30, 2009
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
#include "slu_ddefs.h"
kusano 7d535a
kusano 7d535a
extern int_t mc64id_(int_t*);
kusano 7d535a
extern int_t mc64ad_(int_t*, int_t*, int_t*, int_t [], int_t [], double [],
kusano 7d535a
		    int_t*, int_t [], int_t*, int_t[], int_t*, double [],
kusano 7d535a
		    int_t [], int_t []);
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * Purpose
kusano 7d535a
 * =======
kusano 7d535a
 *
kusano 7d535a
 *   DLDPERM finds a row permutation so that the matrix has large
kusano 7d535a
 *   entries on the diagonal.
kusano 7d535a
 *
kusano 7d535a
 * Arguments
kusano 7d535a
 * =========
kusano 7d535a
 *
kusano 7d535a
 * job    (input) int
kusano 7d535a
 *        Control the action. Possible values for JOB are:
kusano 7d535a
 *        = 1 : Compute a row permutation of the matrix so that the
kusano 7d535a
 *              permuted matrix has as many entries on its diagonal as
kusano 7d535a
 *              possible. The values on the diagonal are of arbitrary size.
kusano 7d535a
 *              HSL subroutine MC21A/AD is used for this.
kusano 7d535a
 *        = 2 : Compute a row permutation of the matrix so that the smallest 
kusano 7d535a
 *              value on the diagonal of the permuted matrix is maximized.
kusano 7d535a
 *        = 3 : Compute a row permutation of the matrix so that the smallest
kusano 7d535a
 *              value on the diagonal of the permuted matrix is maximized.
kusano 7d535a
 *              The algorithm differs from the one used for JOB = 2 and may
kusano 7d535a
 *              have quite a different performance.
kusano 7d535a
 *        = 4 : Compute a row permutation of the matrix so that the sum
kusano 7d535a
 *              of the diagonal entries of the permuted matrix is maximized.
kusano 7d535a
 *        = 5 : Compute a row permutation of the matrix so that the product
kusano 7d535a
 *              of the diagonal entries of the permuted matrix is maximized
kusano 7d535a
 *              and vectors to scale the matrix so that the nonzero diagonal 
kusano 7d535a
 *              entries of the permuted matrix are one in absolute value and 
kusano 7d535a
 *              all the off-diagonal entries are less than or equal to one in 
kusano 7d535a
 *              absolute value.
kusano 7d535a
 *        Restriction: 1 <= JOB <= 5.
kusano 7d535a
 *
kusano 7d535a
 * n      (input) int
kusano 7d535a
 *        The order of the matrix.
kusano 7d535a
 *
kusano 7d535a
 * nnz    (input) int
kusano 7d535a
 *        The number of nonzeros in the matrix.
kusano 7d535a
 *
kusano 7d535a
 * adjncy (input) int*, of size nnz
kusano 7d535a
 *        The adjacency structure of the matrix, which contains the row
kusano 7d535a
 *        indices of the nonzeros.
kusano 7d535a
 *
kusano 7d535a
 * colptr (input) int*, of size n+1
kusano 7d535a
 *        The pointers to the beginning of each column in ADJNCY.
kusano 7d535a
 *
kusano 7d535a
 * nzval  (input) double*, of size nnz
kusano 7d535a
 *        The nonzero values of the matrix. nzval[k] is the value of
kusano 7d535a
 *        the entry corresponding to adjncy[k].
kusano 7d535a
 *        It is not used if job = 1.
kusano 7d535a
 *
kusano 7d535a
 * perm   (output) int*, of size n
kusano 7d535a
 *        The permutation vector. perm[i] = j means row i in the
kusano 7d535a
 *        original matrix is in row j of the permuted matrix.
kusano 7d535a
 *
kusano 7d535a
 * u      (output) double*, of size n
kusano 7d535a
 *        If job = 5, the natural logarithms of the row scaling factors. 
kusano 7d535a
 *
kusano 7d535a
 * v      (output) double*, of size n
kusano 7d535a
 *        If job = 5, the natural logarithms of the column scaling factors. 
kusano 7d535a
 *        The scaled matrix B has entries b_ij = a_ij * exp(u_i + v_j).
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
int
kusano 7d535a
dldperm(int_t job, int_t n, int_t nnz, int_t colptr[], int_t adjncy[],
kusano 7d535a
	double nzval[], int_t *perm, double u[], double v[])
kusano 7d535a
{ 
kusano 7d535a
    int_t i, liw, ldw, num;
kusano 7d535a
    int_t *iw, icntl[10], info[10];
kusano 7d535a
    double *dw;
kusano 7d535a
kusano 7d535a
#if ( DEBUGlevel>=1 )
kusano 7d535a
    CHECK_MALLOC(0, "Enter dldperm()");
kusano 7d535a
#endif
kusano 7d535a
    liw = 5*n;
kusano 7d535a
    if ( job == 3 ) liw = 10*n + nnz;
kusano 7d535a
    if ( !(iw = intMalloc(liw)) ) ABORT("Malloc fails for iw[]");
kusano 7d535a
    ldw = 3*n + nnz;
kusano 7d535a
    if ( !(dw = (double*) SUPERLU_MALLOC(ldw * sizeof(double))) )
kusano 7d535a
          ABORT("Malloc fails for dw[]");
kusano 7d535a
	    
kusano 7d535a
    /* Increment one to get 1-based indexing. */
kusano 7d535a
    for (i = 0; i <= n; ++i) ++colptr[i];
kusano 7d535a
    for (i = 0; i < nnz; ++i) ++adjncy[i];
kusano 7d535a
#if ( DEBUGlevel>=2 )
kusano 7d535a
    printf("LDPERM(): n %d, nnz %d\n", n, nnz);
kusano 7d535a
    slu_PrintInt10("colptr", n+1, colptr);
kusano 7d535a
    slu_PrintInt10("adjncy", nnz, adjncy);
kusano 7d535a
#endif
kusano 7d535a
	
kusano 7d535a
    /* 
kusano 7d535a
     * NOTE:
kusano 7d535a
     * =====
kusano 7d535a
     *
kusano 7d535a
     * MC64AD assumes that column permutation vector is defined as:
kusano 7d535a
     * perm(i) = j means column i of permuted A is in column j of original A.
kusano 7d535a
     *
kusano 7d535a
     * Since a symmetric permutation preserves the diagonal entries. Then
kusano 7d535a
     * by the following relation:
kusano 7d535a
     *     P'(A*P')P = P'A
kusano 7d535a
     * we can apply inverse(perm) to rows of A to get large diagonal entries.
kusano 7d535a
     * But, since 'perm' defined in MC64AD happens to be the reverse of
kusano 7d535a
     * SuperLU's definition of permutation vector, therefore, it is already
kusano 7d535a
     * an inverse for our purpose. We will thus use it directly.
kusano 7d535a
     *
kusano 7d535a
     */
kusano 7d535a
    mc64id_(icntl);
kusano 7d535a
#if 0
kusano 7d535a
    /* Suppress error and warning messages. */
kusano 7d535a
    icntl[0] = -1;
kusano 7d535a
    icntl[1] = -1;
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    mc64ad_(&job, &n, &nnz, colptr, adjncy, nzval, &num, perm,
kusano 7d535a
	    &liw, iw, &ldw, dw, icntl, info);
kusano 7d535a
kusano 7d535a
#if ( DEBUGlevel>=2 )
kusano 7d535a
    slu_PrintInt10("perm", n, perm);
kusano 7d535a
    printf(".. After MC64AD info %d\tsize of matching %d\n", info[0], num);
kusano 7d535a
#endif
kusano 7d535a
    if ( info[0] == 1 ) { /* Structurally singular */
kusano 7d535a
        printf(".. The last %d permutations:\n", n-num);
kusano 7d535a
	slu_PrintInt10("perm", n-num, &perm[num]);
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Restore to 0-based indexing. */
kusano 7d535a
    for (i = 0; i <= n; ++i) --colptr[i];
kusano 7d535a
    for (i = 0; i < nnz; ++i) --adjncy[i];
kusano 7d535a
    for (i = 0; i < n; ++i) --perm[i];
kusano 7d535a
kusano 7d535a
    if ( job == 5 )
kusano 7d535a
        for (i = 0; i < n; ++i) {
kusano 7d535a
	    u[i] = dw[i];
kusano 7d535a
	    v[i] = dw[n+i];
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE(iw);
kusano 7d535a
    SUPERLU_FREE(dw);
kusano 7d535a
kusano 7d535a
#if ( DEBUGlevel>=1 )
kusano 7d535a
    CHECK_MALLOC(0, "Exit dldperm()");
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    return info[0];
kusano 7d535a
}