Actual source code: essl.c

      1: #define PETSCMAT_DLL

      3: /*
      4:         Provides an interface to the IBM RS6000 Essl sparse solver

      6: */
     7:  #include src/mat/impls/aij/seq/aij.h
      8: /* #include <essl.h> This doesn't work!  */

     11: void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*);
     12: void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*);

     15: typedef struct {
     16:   int         n,nz;
     17:   PetscScalar *a;
     18:   int         *ia;
     19:   int         *ja;
     20:   int         lna;
     21:   int         iparm[5];
     22:   PetscReal   rparm[5];
     23:   PetscReal   oparm[5];
     24:   PetscScalar *aux;
     25:   int         naux;

     27:   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
     28:   PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType);
     29:   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
     30:   PetscErrorCode (*MatDestroy)(Mat);
     31:   PetscTruth CleanUpESSL;
     32: } Mat_Essl;

     34: EXTERN PetscErrorCode MatDuplicate_Essl(Mat,MatDuplicateOption,Mat*);

     39: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_Essl_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) {
     41:   Mat            B=*newmat;
     42:   Mat_Essl       *essl=(Mat_Essl*)A->spptr;
     43:
     45:   if (reuse == MAT_INITIAL_MATRIX) {
     46:     MatDuplicate(A,MAT_COPY_VALUES,&B);
     47:   }
     48:   B->ops->duplicate        = essl->MatDuplicate;
     49:   B->ops->assemblyend      = essl->MatAssemblyEnd;
     50:   B->ops->lufactorsymbolic = essl->MatLUFactorSymbolic;
     51:   B->ops->destroy          = essl->MatDestroy;

     53:   /* free the Essl datastructures */
     54:   PetscFree(essl);

     56:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_essl_C","",PETSC_NULL);
     57:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_essl_seqaij_C","",PETSC_NULL);

     59:   PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
     60:   *newmat = B;
     61:   return(0);
     62: }

     67: PetscErrorCode MatDestroy_Essl(Mat A)
     68: {
     70:   Mat_Essl       *essl=(Mat_Essl*)A->spptr;

     73:   if (essl->CleanUpESSL) {
     74:     PetscFree(essl->a);
     75:   }
     76:   MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
     77:   (*A->ops->destroy)(A);
     78:   return(0);
     79: }

     83: PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) {
     84:   Mat_Essl       *essl = (Mat_Essl*)A->spptr;
     85:   PetscScalar    *xx;
     87:   int            m,zero = 0;

     90:   VecGetLocalSize(b,&m);
     91:   VecCopy(b,x);
     92:   VecGetArray(x,&xx);
     93:   dgss(&zero,&A->n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux);
     94:   VecRestoreArray(x,&xx);
     95:   return(0);
     96: }

    100: PetscErrorCode MatLUFactorNumeric_Essl(Mat A,MatFactorInfo *info,Mat *F) {
    101:   Mat_SeqAIJ     *aa=(Mat_SeqAIJ*)(A)->data;
    102:   Mat_Essl       *essl=(Mat_Essl*)(*F)->spptr;
    104:   int            i,one = 1;

    107:   /* copy matrix data into silly ESSL data structure (1-based Frotran style) */
    108:   for (i=0; i<A->m+1; i++) essl->ia[i] = aa->i[i] + 1;
    109:   for (i=0; i<aa->nz; i++) essl->ja[i]  = aa->j[i] + 1;
    110:
    111:   PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));
    112:
    113:   /* set Essl options */
    114:   essl->iparm[0] = 1;
    115:   essl->iparm[1] = 5;
    116:   essl->iparm[2] = 1;
    117:   essl->iparm[3] = 0;
    118:   essl->rparm[0] = 1.e-12;
    119:   essl->rparm[1] = 1.0;
    120:   PetscOptionsGetReal(A->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);

    122:   dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,
    123:                essl->rparm,essl->oparm,essl->aux,&essl->naux);

    125:   (*F)->assembled = PETSC_TRUE;
    126:   return(0);
    127: }

    131: PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) {
    132:   Mat            B;
    133:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
    135:   int            len;
    136:   Mat_Essl       *essl;
    137:   PetscReal      f = 1.0;

    140:   if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
    141:   MatCreate(A->comm,&B);
    142:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n);
    143:   MatSetType(B,A->type_name);
    144:   MatSeqAIJSetPreallocation(B,0,PETSC_NULL);

    146:   B->ops->solve           = MatSolve_Essl;
    147:   B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
    148:   B->factor               = FACTOR_LU;
    149:
    150:   essl = (Mat_Essl *)(B->spptr);

    152:   /* allocate the work arrays required by ESSL */
    153:   f = info->fill;
    154:   essl->nz   = a->nz;
    155:   essl->lna  = (int)a->nz*f;
    156:   essl->naux = 100 + 10*A->m;

    158:   /* since malloc is slow on IBM we try a single malloc */
    159:   len               = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar);
    160:   PetscMalloc(len,&essl->a);
    161:   essl->aux         = essl->a + essl->lna;
    162:   essl->ia          = (int*)(essl->aux + essl->naux);
    163:   essl->ja          = essl->ia + essl->lna;
    164:   essl->CleanUpESSL = PETSC_TRUE;

    166:   PetscLogObjectMemory(B,len+sizeof(Mat_Essl));
    167:   *F = B;
    168:   return(0);
    169: }

    173: PetscErrorCode MatAssemblyEnd_Essl(Mat A,MatAssemblyType mode)
    174: {
    176:   Mat_Essl       *essl=(Mat_Essl*)(A->spptr);

    179:   (*essl->MatAssemblyEnd)(A,mode);

    181:   essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
    182:   A->ops->lufactorsymbolic  = MatLUFactorSymbolic_Essl;
    183:   PetscLogInfo((0,"MatAssemblyEnd_Essl:Using ESSL for LU factorization and solves\n"));
    184:   return(0);
    185: }

    190: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_Essl(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    191: {
    192:   Mat            B=*newmat;
    194:   Mat_Essl       *essl;

    197:   if (reuse == MAT_INITIAL_MATRIX) {
    198:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    199:   }

    201:   PetscNew(Mat_Essl,&essl);
    202:   essl->MatDuplicate        = A->ops->duplicate;
    203:   essl->MatAssemblyEnd      = A->ops->assemblyend;
    204:   essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
    205:   essl->MatDestroy          = A->ops->destroy;
    206:   essl->CleanUpESSL         = PETSC_FALSE;

    208:   B->spptr                  = (void*)essl;
    209:   B->ops->duplicate         = MatDuplicate_Essl;
    210:   B->ops->assemblyend       = MatAssemblyEnd_Essl;
    211:   B->ops->lufactorsymbolic  = MatLUFactorSymbolic_Essl;
    212:   B->ops->destroy           = MatDestroy_Essl;

    214:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C",
    215:                                            "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);
    216:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C",
    217:                                            "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);
    218:   PetscObjectChangeTypeName((PetscObject)B,type);

    220:   *newmat = B;
    221:   return(0);
    222: }

    227: PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M)
    228: {
    230:   Mat_Essl       *lu = (Mat_Essl *)A->spptr;

    233:   (*lu->MatDuplicate)(A,op,M);
    234:   PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));
    235:   return(0);
    236: }

    238: /*MC
    239:   MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices
    240:   via the external package ESSL.

    242:   If ESSL is installed (see the manual for
    243:   instructions on how to declare the existence of external packages),
    244:   a matrix type can be constructed which invokes ESSL solvers.
    245:   After calling MatCreate(...,A), simply call MatSetType(A,MATESSL).
    246:   This matrix type is only supported for double precision real.

    248:   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is
    249:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    250:   the MATSEQAIJ type without data copy.

    252:   Options Database Keys:
    253: . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions()

    255:    Level: beginner

    257: .seealso: PCLU
    258: M*/

    263: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Essl(Mat A)
    264: {

    268:   /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */
    269:   PetscObjectChangeTypeName((PetscObject)A,MATESSL);
    270:   MatSetType(A,MATSEQAIJ);
    271:   MatConvert_SeqAIJ_Essl(A,MATESSL,MAT_REUSE_MATRIX,&A);
    272:   return(0);
    273: }