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, missing; 10 PetscInt reallocs = 0, i, *ai = a->i, *aj = a->j, am = A->rmap->n, *ui; 11 const PetscInt *rip, *riip; 12 PetscInt j; 13 PetscInt d; 14 PetscInt ncols, *cols, *uj; 15 PetscReal fill = info->fill, levels = info->levels; 16 IS iperm; 17 spbas_matrix Pattern_0, Pattern_P; 18 19 PetscFunctionBegin; 20 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); 21 PetscCall(MatMissingDiagonal(A, &missing, &d)); 22 PetscCheck(!missing, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry %" PetscInt_FMT, d); 23 PetscCall(ISIdentity(perm, &perm_identity)); 24 PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm)); 25 26 /* ICC(0) without matrix ordering: simply copies fill pattern */ 27 if (!levels && perm_identity) { 28 PetscCall(PetscMalloc1(am + 1, &ui)); 29 ui[0] = 0; 30 31 for (i = 0; i < am; i++) ui[i + 1] = ui[i] + ai[i + 1] - a->diag[i]; 32 PetscCall(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 PetscCall(ISGetIndices(iperm, &riip)); 41 PetscCall(ISGetIndices(perm, &rip)); 42 43 /* Create spbas_matrix for pattern */ 44 PetscCall(spbas_pattern_only(am, am, ai, aj, &Pattern_0)); 45 46 /* Apply the permutation */ 47 PetscCall(spbas_apply_reordering(&Pattern_0, rip, riip)); 48 49 /* Raise the power */ 50 PetscCall(spbas_power(Pattern_0, (int)levels + 1, &Pattern_P)); 51 PetscCall(spbas_delete(Pattern_0)); 52 53 /* Keep only upper triangle of pattern */ 54 PetscCall(spbas_keep_upper(&Pattern_P)); 55 56 /* Convert to Sparse Row Storage */ 57 PetscCall(spbas_matrix_to_crs(Pattern_P, NULL, &ui, &uj)); 58 PetscCall(spbas_delete(Pattern_P)); 59 } /* end of case: levels>0 || (levels=0 && !perm_identity) */ 60 61 /* put together the new matrix in MATSEQSBAIJ format */ 62 63 b = (Mat_SeqSBAIJ *)(fact)->data; 64 b->singlemalloc = PETSC_FALSE; 65 66 PetscCall(PetscMalloc1(ui[am] + 1, &b->a)); 67 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; 75 76 PetscCall(PetscObjectReference((PetscObject)perm)); 77 PetscCall(PetscObjectReference((PetscObject)perm)); 78 79 b->icol = iperm; 80 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 81 PetscCall(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; 85 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 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 static 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; 110 111 PetscFunctionBegin; 112 /* Reduce memory requirements: erase values of B-matrix */ 113 PetscCall(PetscFree(ba)); 114 /* Compress (maximum) sparseness pattern of B-matrix */ 115 PetscCall(spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS, &Pattern, &mem_reduction)); 116 PetscCall(PetscFree(bi)); 117 PetscCall(PetscFree(bj)); 118 119 PetscCall(PetscInfo(NULL, " compression rate for spbas_compress_pattern %g \n", (double)mem_reduction)); 120 121 /* Make Cholesky decompositions with larger Manteuffel shifts until no more negative diagonals are found. */ 122 PetscCall(ISGetIndices(ip, &rip)); 123 PetscCall(ISGetIndices(iip, &riip)); 124 125 if (info->usedt) droptol = info->dt; 126 127 for (int ierr = NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) { 128 PetscBool success; 129 130 ierr = (int)spbas_incomplete_cholesky(A, rip, riip, Pattern, droptol, shiftnz, &matrix_LT, &success); 131 if (!success) { 132 shiftnz *= 1.5; 133 if (shiftnz < 1e-5) shiftnz = 1e-5; 134 PetscCall(PetscInfo(NULL, "spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%g\n", (double)shiftnz)); 135 } 136 } 137 PetscCall(spbas_delete(Pattern)); 138 139 PetscCall(PetscInfo(NULL, " memory_usage for spbas_incomplete_cholesky %g bytes per row\n", (double)(PetscReal)(spbas_memory_requirement(matrix_LT) / (PetscReal)mbs))); 140 141 PetscCall(ISRestoreIndices(ip, &rip)); 142 PetscCall(ISRestoreIndices(iip, &riip)); 143 144 /* Convert spbas_matrix to compressed row storage */ 145 PetscCall(spbas_transpose(matrix_LT, &matrix_L)); 146 PetscCall(spbas_delete(matrix_LT)); 147 PetscCall(spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj)); 148 b->i = bi; 149 b->j = bj; 150 b->a = ba; 151 PetscCall(spbas_delete(matrix_L)); 152 153 /* Set the appropriate solution functions */ 154 PetscCall(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 } 166 167 C->assembled = PETSC_TRUE; 168 C->preallocated = PETSC_TRUE; 169 170 PetscCall(PetscLogFlops(C->rmap->n)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 174 static PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A, MatSolverType *type) 175 { 176 PetscFunctionBegin; 177 *type = MATSOLVERBAS; 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat A, MatFactorType ftype, Mat *B) 182 { 183 PetscInt n = A->rmap->n; 184 185 PetscFunctionBegin; 186 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 187 PetscCall(MatSetSizes(*B, n, n, n, n)); 188 if (ftype == MAT_FACTOR_ICC) { 189 PetscCall(MatSetType(*B, MATSEQSBAIJ)); 190 PetscCall(MatSeqSBAIJSetPreallocation(*B, 1, MAT_SKIP_ALLOCATION, NULL)); 191 192 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqAIJ_Bas; 193 (*B)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_Bas; 194 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_bas)); 195 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU])); 196 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 197 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 198 (*B)->factortype = ftype; 199 200 PetscCall(PetscFree((*B)->solvertype)); 201 PetscCall(PetscStrallocpy(MATSOLVERBAS, &(*B)->solvertype)); 202 (*B)->canuseordering = PETSC_TRUE; 203 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC])); 204 PetscFunctionReturn(PETSC_SUCCESS); 205 } 206