1 #include <../src/mat/impls/aij/seq/aij.h>
2 #include <../src/mat/impls/sbaij/seq/sbaij.h>
3 #include <../src/mat/impls/aij/seq/bas/spbas.h>
4
MatICCFactorSymbolic_SeqAIJ_Bas(Mat fact,Mat A,IS perm,const MatFactorInfo * info)5 static PetscErrorCode MatICCFactorSymbolic_SeqAIJ_Bas(Mat fact, Mat A, IS perm, const MatFactorInfo *info)
6 {
7 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
8 Mat_SeqSBAIJ *b;
9 PetscBool perm_identity, diagDense;
10 PetscInt reallocs = 0, i, *ai = a->i, *aj = a->j, am = A->rmap->n, *ui;
11 const PetscInt *rip, *riip, *adiag;
12 PetscInt j;
13 PetscInt ncols, *cols, *uj;
14 PetscReal fill = info->fill, levels = info->levels;
15 IS iperm;
16 spbas_matrix Pattern_0, Pattern_P;
17
18 PetscFunctionBegin;
19 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, A->rmap->n, A->cmap->n);
20 PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, &diagDense));
21 PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entries");
22 PetscCall(ISIdentity(perm, &perm_identity));
23 PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
24
25 /* ICC(0) without matrix ordering: simply copies fill pattern */
26 if (!levels && perm_identity) {
27 PetscCall(PetscMalloc1(am + 1, &ui));
28 ui[0] = 0;
29
30 for (i = 0; i < am; i++) ui[i + 1] = ui[i] + ai[i + 1] - adiag[i];
31 PetscCall(PetscMalloc1(ui[am] + 1, &uj));
32 cols = uj;
33 for (i = 0; i < am; i++) {
34 aj = a->j + adiag[i];
35 ncols = ui[i + 1] - ui[i];
36 for (j = 0; j < ncols; j++) *cols++ = *aj++;
37 }
38 } else { /* case: levels>0 || (levels=0 && !perm_identity) */
39 PetscCall(ISGetIndices(iperm, &riip));
40 PetscCall(ISGetIndices(perm, &rip));
41
42 /* Create spbas_matrix for pattern */
43 PetscCall(spbas_pattern_only(am, am, ai, aj, &Pattern_0));
44
45 /* Apply the permutation */
46 PetscCall(spbas_apply_reordering(&Pattern_0, rip, riip));
47
48 /* Raise the power */
49 PetscCall(spbas_power(Pattern_0, (int)levels + 1, &Pattern_P));
50 PetscCall(spbas_delete(Pattern_0));
51
52 /* Keep only upper triangle of pattern */
53 PetscCall(spbas_keep_upper(&Pattern_P));
54
55 /* Convert to Sparse Row Storage */
56 PetscCall(spbas_matrix_to_crs(Pattern_P, NULL, &ui, &uj));
57 PetscCall(spbas_delete(Pattern_P));
58 } /* end of case: levels>0 || (levels=0 && !perm_identity) */
59
60 /* put together the new matrix in MATSEQSBAIJ format */
61
62 b = (Mat_SeqSBAIJ *)fact->data;
63 PetscCall(PetscMalloc1(ui[am], &b->a));
64
65 b->j = uj;
66 b->i = ui;
67 b->diag = NULL;
68 b->ilen = NULL;
69 b->imax = NULL;
70 b->row = perm;
71 b->col = perm;
72
73 PetscCall(PetscObjectReference((PetscObject)perm));
74 PetscCall(PetscObjectReference((PetscObject)perm));
75
76 b->icol = iperm;
77 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
78 PetscCall(PetscMalloc1(am, &b->solve_work));
79 b->maxnz = b->nz = ui[am];
80 b->free_a = PETSC_TRUE;
81 b->free_ij = PETSC_TRUE;
82
83 fact->info.factor_mallocs = reallocs;
84 fact->info.fill_ratio_given = fill;
85 if (ai[am] != 0) {
86 fact->info.fill_ratio_needed = (PetscReal)ui[am] / (PetscReal)ai[am];
87 } else {
88 fact->info.fill_ratio_needed = 0.0;
89 }
90 /* fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace; */
91 PetscFunctionReturn(PETSC_SUCCESS);
92 }
93
MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B,Mat A,const MatFactorInfo * info)94 static PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B, Mat A, const MatFactorInfo *info)
95 {
96 Mat C = B;
97 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data;
98 IS ip = b->row, iip = b->icol;
99 const PetscInt *rip, *riip;
100 PetscInt mbs = A->rmap->n, *bi = b->i, *bj = b->j;
101 MatScalar *ba = b->a;
102 PetscReal shiftnz = info->shiftamount;
103 PetscReal droptol = -1;
104 PetscBool perm_identity;
105 spbas_matrix Pattern, matrix_L, matrix_LT;
106 PetscReal mem_reduction;
107
108 PetscFunctionBegin;
109 /* Reduce memory requirements: erase values of B-matrix */
110 PetscCall(PetscFree(ba));
111 /* Compress (maximum) sparseness pattern of B-matrix */
112 PetscCall(spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS, &Pattern, &mem_reduction));
113 PetscCall(PetscFree(bi));
114 PetscCall(PetscFree(bj));
115
116 PetscCall(PetscInfo(NULL, " compression rate for spbas_compress_pattern %g \n", (double)mem_reduction));
117
118 /* Make Cholesky decompositions with larger Manteuffel shifts until no more negative diagonals are found. */
119 PetscCall(ISGetIndices(ip, &rip));
120 PetscCall(ISGetIndices(iip, &riip));
121
122 if (info->usedt) droptol = info->dt;
123
124 for (int ierr = NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) {
125 PetscBool success;
126
127 ierr = (int)spbas_incomplete_cholesky(A, rip, riip, Pattern, droptol, shiftnz, &matrix_LT, &success);
128 if (!success) {
129 shiftnz *= 1.5;
130 if (shiftnz < 1e-5) shiftnz = 1e-5;
131 PetscCall(PetscInfo(NULL, "spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%g\n", (double)shiftnz));
132 }
133 }
134 PetscCall(spbas_delete(Pattern));
135
136 PetscCall(PetscInfo(NULL, " memory_usage for spbas_incomplete_cholesky %g bytes per row\n", (double)(spbas_memory_requirement(matrix_LT) / (PetscReal)mbs)));
137
138 PetscCall(ISRestoreIndices(ip, &rip));
139 PetscCall(ISRestoreIndices(iip, &riip));
140
141 /* Convert spbas_matrix to compressed row storage */
142 PetscCall(spbas_transpose(matrix_LT, &matrix_L));
143 PetscCall(spbas_delete(matrix_LT));
144 PetscCall(spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj));
145 b->i = bi;
146 b->j = bj;
147 b->a = ba;
148 PetscCall(spbas_delete(matrix_L));
149
150 /* Set the appropriate solution functions */
151 PetscCall(ISIdentity(ip, &perm_identity));
152 if (perm_identity) {
153 B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
154 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
155 B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
156 B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace;
157 } else {
158 B->ops->solve = MatSolve_SeqSBAIJ_1_inplace;
159 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace;
160 B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace;
161 B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace;
162 }
163
164 C->assembled = PETSC_TRUE;
165 C->preallocated = PETSC_TRUE;
166
167 PetscCall(PetscLogFlops(C->rmap->n));
168 PetscFunctionReturn(PETSC_SUCCESS);
169 }
170
MatFactorGetSolverType_seqaij_bas(Mat A,MatSolverType * type)171 static PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A, MatSolverType *type)
172 {
173 PetscFunctionBegin;
174 *type = MATSOLVERBAS;
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
MatGetFactor_seqaij_bas(Mat A,MatFactorType ftype,Mat * B)178 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat A, MatFactorType ftype, Mat *B)
179 {
180 PetscInt n = A->rmap->n;
181
182 PetscFunctionBegin;
183 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
184 PetscCall(MatSetSizes(*B, n, n, n, n));
185 PetscCheck(ftype == MAT_FACTOR_ICC, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
186 PetscCall(MatSetType(*B, MATSEQSBAIJ));
187 PetscCall(MatSeqSBAIJSetPreallocation(*B, 1, MAT_SKIP_ALLOCATION, NULL));
188
189 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqAIJ_Bas;
190 (*B)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_Bas;
191 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_bas));
192 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU]));
193 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]));
194 (*B)->factortype = ftype;
195
196 PetscCall(PetscFree((*B)->solvertype));
197 PetscCall(PetscStrallocpy(MATSOLVERBAS, &(*B)->solvertype));
198 (*B)->canuseordering = PETSC_TRUE;
199 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC]));
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202