Actual source code: basfactor.c


  2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
  4: #include <../src/mat/impls/aij/seq/bas/spbas.h>

  6: PetscErrorCode MatICCFactorSymbolic_SeqAIJ_Bas(Mat fact,Mat A,IS perm,const MatFactorInfo *info)
  7: {
  8:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
  9:   Mat_SeqSBAIJ   *b;
 10:   PetscBool      perm_identity,missing;
 11:   PetscInt       reallocs=0,i,*ai=a->i,*aj=a->j,am=A->rmap->n,*ui;
 12:   const PetscInt *rip,*riip;
 13:   PetscInt       j;
 14:   PetscInt       d;
 15:   PetscInt       ncols,*cols,*uj;
 16:   PetscReal      fill=info->fill,levels=info->levels;
 17:   IS             iperm;
 18:   spbas_matrix   Pattern_0, Pattern_P;

 21:   MatMissingDiagonal(A,&missing,&d);
 23:   ISIdentity(perm,&perm_identity);
 24:   ISInvertPermutation(perm,PETSC_DECIDE,&iperm);

 26:   /* ICC(0) without matrix ordering: simply copies fill pattern */
 27:   if (!levels && perm_identity) {
 28:     PetscMalloc1(am+1,&ui);
 29:     ui[0] = 0;

 31:     for (i=0; i<am; i++) {
 32:       ui[i+1] = ui[i] + ai[i+1] - a->diag[i];
 33:     }
 34:     PetscMalloc1(ui[am]+1,&uj);
 35:     cols = uj;
 36:     for (i=0; i<am; i++) {
 37:       aj    = a->j + a->diag[i];
 38:       ncols = ui[i+1] - ui[i];
 39:       for (j=0; j<ncols; j++) *cols++ = *aj++;
 40:     }
 41:   } else { /* case: levels>0 || (levels=0 && !perm_identity) */
 42:     ISGetIndices(iperm,&riip);
 43:     ISGetIndices(perm,&rip);

 45:     /* Create spbas_matrix for pattern */
 46:     spbas_pattern_only(am, am, ai, aj, &Pattern_0);

 48:     /* Apply the permutation */
 49:     spbas_apply_reordering(&Pattern_0, rip, riip);

 51:     /* Raise the power */
 52:     spbas_power(Pattern_0, (int) levels+1, &Pattern_P);
 53:     spbas_delete(Pattern_0);

 55:     /* Keep only upper triangle of pattern */
 56:     spbas_keep_upper(&Pattern_P);

 58:     /* Convert to Sparse Row Storage  */
 59:     spbas_matrix_to_crs(Pattern_P, NULL, &ui, &uj);
 60:     spbas_delete(Pattern_P);
 61:   } /* end of case: levels>0 || (levels=0 && !perm_identity) */

 63:   /* put together the new matrix in MATSEQSBAIJ format */

 65:   b               = (Mat_SeqSBAIJ*)(fact)->data;
 66:   b->singlemalloc = PETSC_FALSE;

 68:   PetscMalloc1(ui[am]+1,&b->a);

 70:   b->j    = uj;
 71:   b->i    = ui;
 72:   b->diag = NULL;
 73:   b->ilen = NULL;
 74:   b->imax = NULL;
 75:   b->row  = perm;
 76:   b->col  = perm;

 78:   PetscObjectReference((PetscObject)perm);
 79:   PetscObjectReference((PetscObject)perm);

 81:   b->icol          = iperm;
 82:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
 83:   PetscMalloc1(am+1,&b->solve_work);
 84:   PetscLogObjectMemory((PetscObject)(fact),(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));
 85:   b->maxnz         = b->nz = ui[am];
 86:   b->free_a        = PETSC_TRUE;
 87:   b->free_ij       = PETSC_TRUE;

 89:   (fact)->info.factor_mallocs   = reallocs;
 90:   (fact)->info.fill_ratio_given = fill;
 91:   if (ai[am] != 0) {
 92:     (fact)->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]);
 93:   } else {
 94:     (fact)->info.fill_ratio_needed = 0.0;
 95:   }
 96:   /*  (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace; */
 97:   return 0;
 98: }

100: PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B,Mat A,const MatFactorInfo *info)
101: {
102:   Mat            C = B;
103:   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
104:   IS             ip=b->row,iip = b->icol;
105:   const PetscInt *rip,*riip;
106:   PetscInt       mbs=A->rmap->n,*bi=b->i,*bj=b->j;
107:   MatScalar      *ba     = b->a;
108:   PetscReal      shiftnz = info->shiftamount;
109:   PetscReal      droptol = -1;
110:   PetscBool      perm_identity;
111:   spbas_matrix   Pattern, matrix_L,matrix_LT;
112:   PetscReal      mem_reduction;

114:   /* Reduce memory requirements:   erase values of B-matrix */
115:   PetscFree(ba);
116:   /*   Compress (maximum) sparseness pattern of B-matrix */
117:   spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS,&Pattern, &mem_reduction);
118:   PetscFree(bi);
119:   PetscFree(bj);

121:   PetscInfo(NULL,"    compression rate for spbas_compress_pattern %g \n",(double)mem_reduction);

123:   /* Make Cholesky decompositions with larger Manteuffel shifts until no more    negative diagonals are found. */
124:   ISGetIndices(ip,&rip);
125:   ISGetIndices(iip,&riip);

127:   if (info->usedt) droptol = info->dt;

129:   for (PetscErrorCode NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) {
130:     PetscBool success;

132:     spbas_incomplete_cholesky(A, rip, riip, Pattern, droptol, shiftnz,&matrix_LT,&success);
133:     if (!success) {
134:       shiftnz *= 1.5;
135:       if (shiftnz < 1e-5) shiftnz=1e-5;
136:       PetscInfo(NULL,"spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%g\n",(double)shiftnz);
137:     }
138:   }
139:   spbas_delete(Pattern);

141:   PetscInfo(NULL,"    memory_usage for  spbas_incomplete_cholesky  %g bytes per row\n", (double)(PetscReal) (spbas_memory_requirement(matrix_LT)/ (PetscReal) mbs));

143:   ISRestoreIndices(ip,&rip);
144:   ISRestoreIndices(iip,&riip);

146:   /* Convert spbas_matrix to compressed row storage */
147:   spbas_transpose(matrix_LT, &matrix_L);
148:   spbas_delete(matrix_LT);
149:   spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj);
150:   b->i =bi; b->j=bj; b->a=ba;
151:   spbas_delete(matrix_L);

153:   /* Set the appropriate solution functions */
154:   ISIdentity(ip,&perm_identity);
155:   if (perm_identity) {
156:     (B)->ops->solve          = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
157:     (B)->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
158:     (B)->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
159:     (B)->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
160:   } else {
161:     (B)->ops->solve          = MatSolve_SeqSBAIJ_1_inplace;
162:     (B)->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
163:     (B)->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_1_inplace;
164:     (B)->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_1_inplace;
165:   }

167:   C->assembled    = PETSC_TRUE;
168:   C->preallocated = PETSC_TRUE;

170:   PetscLogFlops(C->rmap->n);
171:   return 0;
172: }

174: PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A,MatSolverType *type)
175: {
176:   *type = MATSOLVERBAS;
177:   return 0;
178: }

180: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat A,MatFactorType ftype,Mat *B)
181: {
182:   PetscInt       n = A->rmap->n;

184:   MatCreate(PetscObjectComm((PetscObject)A),B);
185:   MatSetSizes(*B,n,n,n,n);
186:   if (ftype == MAT_FACTOR_ICC) {
187:     MatSetType(*B,MATSEQSBAIJ);
188:     MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,NULL);

190:     (*B)->ops->iccfactorsymbolic     = MatICCFactorSymbolic_SeqAIJ_Bas;
191:     (*B)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_Bas;
192:     PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_bas);
193:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_LU]);
194:     PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);
195:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
196:   (*B)->factortype = ftype;

198:   PetscFree((*B)->solvertype);
199:   PetscStrallocpy(MATSOLVERBAS,&(*B)->solvertype);
200:   (*B)->canuseordering = PETSC_TRUE;
201:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);
202:   return 0;
203: }