superlu.c

Go to the documentation of this file.

/* Wrapper for SuperLU solver. Based on fortran wrapper supplied
 * with SuperLU version 3.0. Given that, it seems appropriate
 * to retain this copyright notice:
 *
 * -- SuperLU routine (version 3.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * October 15, 2003
 *
 */

#ifdef USING_OOMPH_SUPERLU
#include "../../external_src/oomph_superlu_4.3/slu_ddefs.h"
#else
#include "slu_ddefs.h"
#endif
#include "math.h"

/* ================================================= */
/* Pointer to the LU factors*/
/* ================================================= */
typedef void* fptr;  


/* ================================================= */
/* Struct for the lu factors  */
/* ================================================= */
typedef struct {
    SuperMatrix *L;
    SuperMatrix *U;
    int *perm_c;
    int *perm_r;
} factors_t;




/* =========================================================================
 * Wrapper to superlu solver:
 *
 * op_flag    = int specifies the operation:
 *                1, performs LU decomposition for the first time
 *                2, performs triangular solve
 *                3, free all the storage in the end
 * n          = dimension of matrix
 * nnz        = # of nonzero entries
 * nrhs       = # of RHSs
 * values     = double array containing the nonzero entries 
 * rowind     = int array containing the row indices of the entries 
 * colptr     = int array containing the column start
 * b          = double array containing the rhs vector (gets overwritten
 *              with solution)
 * ldb        = leading dimension of b
 * transpose  =  0/1 if matrix is transposed/not transposed
 * doc        = 0/1 for full doc/no full doc 
 * info       = info flag from superlu
 * f_factors  = pointer to LU factors. (If op_flag == 1, it is an output 
 *              and contains the pointer pointing to the structure of 
 *              the factored matrices. Otherwise, it it an input.
 * Returns the SIGN of the determinant of the matrix 
 * =========================================================================
 */
int superlu(int *op_flag, int *n, int *nnz, int *nrhs,
             double *values, int *rowind, int *colptr,
             double *b, int *ldb, int *transpose, int *doc,
             fptr *f_factors, int *info)
 
{
 
    SuperMatrix A, AC, B;
    SuperMatrix *L, *U;
    int *perm_r; /* row permutations from partial pivoting */
    int *perm_c; /* column permutation vector */
    int *etree;  /* column elimination tree */
    SCformat *Lstore;
    NCformat *Ustore;
    int      i, j, panel_size, permc_spec, relax;
    trans_t  trans;
    //double   drop_tol = 0.0; //No longer need SuperLU 4.3
    mem_usage_t   mem_usage;
    superlu_options_t options;
    SuperLUStat_t stat;
    factors_t *LUfactors;

    double *Lval;
    double *diagU, *dblock;
    int_t fsupc, nsupr, nsupc, luptr;
    int_t i2, k2, nsupers; 
    int signature=1;
    int sign = 0;

/*   Do we need to transpose? */
    if (*transpose==0)
     {
      trans = NOTRANS;
     }
    else
     {
      trans = TRANS;
     }

    if ( *op_flag == 1 ) { /* LU decomposition */

        /* Set the default input options. */
        set_default_options(&options);

        /* Initialize the statistics variables. */
        StatInit(&stat);

        dCreate_CompCol_Matrix(&A, *n, *n, *nnz, values, rowind, colptr,
                               SLU_NC, SLU_D, SLU_GE);
        L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
        U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
        if ( !(perm_r = intMalloc(*n)) ) ABORT("Malloc fails for perm_r[].");
        if ( !(perm_c = intMalloc(*n)) ) ABORT("Malloc fails for perm_c[].");
        if ( !(etree = intMalloc(*n)) ) ABORT("Malloc fails for etree[].");

        /*
         * Get column permutation vector perm_c[], according to permc_spec:
         *   permc_spec = 0: natural ordering 
         *   permc_spec = 1: minimum degree on structure of A'*A
         *   permc_spec = 2: minimum degree on structure of A'+A
         *   permc_spec = 3: approximate minimum degree for unsymmetric matrices
         */     
        permc_spec = options.ColPerm;
        get_perm_c(permc_spec, &A, perm_c);
        
        sp_preorder(&options, &A, perm_c, etree, &AC);

        panel_size = sp_ienv(1);
        relax = sp_ienv(2);

        dgstrf(&options, &AC, /*drop_tol,*/ relax, panel_size, 
               etree, NULL, 0, perm_c, perm_r, L, U, &stat, info);

        if ( *info == 0 ) {
            Lstore = (SCformat *) L->Store;
            Ustore = (NCformat *) U->Store;
            if (*doc!=0)
             {
              printf(" No of nonzeros in factor L = %d\n", Lstore->nnz);
              printf(" No of nonzeros in factor U = %d\n", Ustore->nnz);
              printf(" No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
              dQuerySpace(L, U, &mem_usage);
              printf(" L\\U MB %.3f\ttotal MB needed %.3f\n",
                     //\texpansions %d\n",
                     mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
                     //mem_usage.expansions);
             }
        }
        else {
            printf("dgstrf() error returns INFO= %d\n", *info);
            if ( *info < 0) {
                printf("Argument number %d had an illegal value", *info);
            }
            else if ( *info <= *n ) { /* factorization completes */
                printf("U(%i,%i) is exactly zero so U is exactly singular.",
                       *n, *n);
               dQuerySpace(L, U, &mem_usage);
               printf(" L\\U MB %.3f\ttotal MB needed %.3f\n",
                      //\texpansions %d\n\n",
                      mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
               //mem_usage.expansions);
            }
        }
        
        /* Save the LU factors in the factors handle */
        LUfactors = (factors_t*) SUPERLU_MALLOC(sizeof(factors_t));
        LUfactors->L = L;
        LUfactors->U = U;
        LUfactors->perm_c = perm_c;
        LUfactors->perm_r = perm_r;
        *f_factors = (fptr) LUfactors;

        //Work out and print the sign of the determinant
        //This code is hacked from supraLU by  Alex Pletzer
        //and Doug McCune (NTCC) (http://w3.pppl.gov/ntcc/SupraLu/)
        Lstore = L->Store;
        Lval = Lstore->nzval;
        nsupers = Lstore->nsuper + 1;

        //Get the diagonal entries of the U matrix
        //Allocate store for the entries
        if ( !(diagU = SUPERLU_MALLOC( *n * sizeof(SuperMatrix) )) )
         ABORT("Malloc fails for diagU[].");
        //Loop over the number of super diagonal terms(?)
        for(k2=0; k2< nsupers; k2++)
         {
          fsupc = L_FST_SUPC(k2);
          nsupc = L_FST_SUPC(k2+1) - fsupc;
          nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
          luptr = L_NZ_START(fsupc);

          dblock = &diagU[fsupc];
          for(i2 = 0; i2 < nsupc; i2++)
           {
            dblock[i2] = Lval[luptr];
            luptr += nsupr + 1;
           }
         }
                  
        //Now multiply all the diagonal terms together to get the determinant
        //Note that we need to use the mantissa, exponent formulation to
        //avoid underflow errors
        double determinant_mantissa=1.0;
        int determinant_exponent = 0, iexp;
        for(i=0; i<*n; i++)
         {
          determinant_mantissa *= frexp(diagU[i], &iexp);
          determinant_exponent += iexp;
          /* normalise*/
          determinant_mantissa = frexp(determinant_mantissa,&iexp);
          determinant_exponent += iexp;
                    
          /*Now worry about the permutations 
            (this is done in a stupid, but not too inefficient way)*/
          for(j=i;j<*n;j++)
           {
            if(perm_r[j] < perm_r[i]) {signature *= -1;}
            if(perm_c[j] < perm_c[i]) {signature *= -1;}
           }
         }
        
        //Find the sign of the determinant
        if(determinant_mantissa > 0.0) {sign = 1;}
        if(determinant_mantissa < 0.0) {sign = -1;}
        
        //Multiply the sign by the signature
        sign *= signature;
 
        /* Free un-wanted storage */
        SUPERLU_FREE(diagU);
        SUPERLU_FREE(etree);
        Destroy_SuperMatrix_Store(&A);
        Destroy_CompCol_Permuted(&AC);
        StatFree(&stat);

        //Return the sign of the determinant
        return sign;
        
    } else if ( *op_flag == 2 ) { /* Triangular solve */
        /* Initialize the statistics variables. */
        StatInit(&stat);

        /* Extract the LU factors in the factors handle */
        LUfactors = (factors_t*) *f_factors;
        L = LUfactors->L;
        U = LUfactors->U;
        perm_c = LUfactors->perm_c;
        perm_r = LUfactors->perm_r;

        dCreate_Dense_Matrix(&B, *n, *nrhs, b, *ldb, SLU_DN, SLU_D, SLU_GE);

        /* Solve the system A*X=B, overwriting B with X. */
        dgstrs (trans, L, U, perm_c, perm_r, &B, &stat, info);

        Destroy_SuperMatrix_Store(&B);
        StatFree(&stat);

        //Return zero
        return 0;

    } else if ( *op_flag == 3 ) { /* Free storage */
        /* Free the LU factors in the factors handle */
        LUfactors = (factors_t*) *f_factors;
        SUPERLU_FREE (LUfactors->perm_r);
        SUPERLU_FREE (LUfactors->perm_c);
        Destroy_SuperNode_Matrix(LUfactors->L);
        Destroy_CompCol_Matrix(LUfactors->U);
        SUPERLU_FREE (LUfactors->L);
        SUPERLU_FREE (LUfactors->U);
        SUPERLU_FREE (LUfactors);
        return 0;
    } else {
        fprintf(stderr,"Invalid op_flag=%d passed to c_cpp_dgssv()\n",*op_flag);
        exit(-1);
    }
}