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 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 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 171 static PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A, MatSolverType *type) 172 { 173 PetscFunctionBegin; 174 *type = MATSOLVERBAS; 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 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