1 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> 5 6 PetscErrorCode MatICCFactorSymbolic_SeqAIJ_Bas(Mat fact, Mat A, IS perm, const MatFactorInfo *info) { 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(0); 95 } 96 97 PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B, Mat A, const MatFactorInfo *info) { 98 Mat C = B; 99 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data; 100 IS ip = b->row, iip = b->icol; 101 const PetscInt *rip, *riip; 102 PetscInt mbs = A->rmap->n, *bi = b->i, *bj = b->j; 103 MatScalar *ba = b->a; 104 PetscReal shiftnz = info->shiftamount; 105 PetscReal droptol = -1; 106 PetscBool perm_identity; 107 spbas_matrix Pattern, matrix_L, matrix_LT; 108 PetscReal mem_reduction; 109 110 PetscFunctionBegin; 111 /* Reduce memory requirements: erase values of B-matrix */ 112 PetscCall(PetscFree(ba)); 113 /* Compress (maximum) sparseness pattern of B-matrix */ 114 PetscCall(spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS, &Pattern, &mem_reduction)); 115 PetscCall(PetscFree(bi)); 116 PetscCall(PetscFree(bj)); 117 118 PetscCall(PetscInfo(NULL, " compression rate for spbas_compress_pattern %g \n", (double)mem_reduction)); 119 120 /* Make Cholesky decompositions with larger Manteuffel shifts until no more negative diagonals are found. */ 121 PetscCall(ISGetIndices(ip, &rip)); 122 PetscCall(ISGetIndices(iip, &riip)); 123 124 if (info->usedt) droptol = info->dt; 125 126 for (PetscErrorCode ierr = NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) { 127 PetscBool success; 128 129 ierr = 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 PetscCall(PetscInfo(NULL, "spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%g\n", (double)shiftnz)); 134 } 135 } 136 PetscCall(spbas_delete(Pattern)); 137 138 PetscCall(PetscInfo(NULL, " memory_usage for spbas_incomplete_cholesky %g bytes per row\n", (double)(PetscReal)(spbas_memory_requirement(matrix_LT) / (PetscReal)mbs))); 139 140 PetscCall(ISRestoreIndices(ip, &rip)); 141 PetscCall(ISRestoreIndices(iip, &riip)); 142 143 /* Convert spbas_matrix to compressed row storage */ 144 PetscCall(spbas_transpose(matrix_LT, &matrix_L)); 145 PetscCall(spbas_delete(matrix_LT)); 146 PetscCall(spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj)); 147 b->i = bi; 148 b->j = bj; 149 b->a = ba; 150 PetscCall(spbas_delete(matrix_L)); 151 152 /* Set the appropriate solution functions */ 153 PetscCall(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 } 165 166 C->assembled = PETSC_TRUE; 167 C->preallocated = PETSC_TRUE; 168 169 PetscCall(PetscLogFlops(C->rmap->n)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A, MatSolverType *type) { 174 PetscFunctionBegin; 175 *type = MATSOLVERBAS; 176 PetscFunctionReturn(0); 177 } 178 179 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat A, MatFactorType ftype, Mat *B) { 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 if (ftype == MAT_FACTOR_ICC) { 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 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 195 (*B)->factortype = ftype; 196 197 PetscCall(PetscFree((*B)->solvertype)); 198 PetscCall(PetscStrallocpy(MATSOLVERBAS, &(*B)->solvertype)); 199 (*B)->canuseordering = PETSC_TRUE; 200 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC])); 201 PetscFunctionReturn(0); 202 } 203