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

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

 46:     /* Apply the permutation */
 47:     spbas_apply_reordering(&Pattern_0, rip, riip);

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

 53:     /* Keep only upper triangle of pattern */
 54:     spbas_keep_upper(&Pattern_P);

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

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

 63:   b               = (Mat_SeqSBAIJ *)(fact)->data;
 64:   b->singlemalloc = PETSC_FALSE;

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

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

 76:   PetscObjectReference((PetscObject)perm);
 77:   PetscObjectReference((PetscObject)perm);

 79:   b->icol          = iperm;
 80:   b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
 81:   PetscMalloc1(am + 1, &b->solve_work);
 82:   b->maxnz = b->nz = ui[am];
 83:   b->free_a        = PETSC_TRUE;
 84:   b->free_ij       = PETSC_TRUE;

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

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

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

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

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

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

126:   for (PetscErrorCode NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) {
127:     PetscBool success;

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

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

140:   ISRestoreIndices(ip, &rip);
141:   ISRestoreIndices(iip, &riip);

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

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

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

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

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

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

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

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

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