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