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: }