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