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 PetscCall(PetscMalloc1(ui[am], &b->a)); 65 66 b->j = uj; 67 b->i = ui; 68 b->diag = NULL; 69 b->ilen = NULL; 70 b->imax = NULL; 71 b->row = perm; 72 b->col = perm; 73 74 PetscCall(PetscObjectReference((PetscObject)perm)); 75 PetscCall(PetscObjectReference((PetscObject)perm)); 76 77 b->icol = iperm; 78 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 79 PetscCall(PetscMalloc1(am, &b->solve_work)); 80 b->maxnz = b->nz = ui[am]; 81 b->free_a = PETSC_TRUE; 82 b->free_ij = PETSC_TRUE; 83 84 fact->info.factor_mallocs = reallocs; 85 fact->info.fill_ratio_given = fill; 86 if (ai[am] != 0) { 87 fact->info.fill_ratio_needed = (PetscReal)ui[am] / (PetscReal)ai[am]; 88 } else { 89 fact->info.fill_ratio_needed = 0.0; 90 } 91 /* fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace; */ 92 PetscFunctionReturn(PETSC_SUCCESS); 93 } 94 95 static PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B, Mat A, const MatFactorInfo *info) 96 { 97 Mat C = B; 98 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)C->data; 99 IS ip = b->row, iip = b->icol; 100 const PetscInt *rip, *riip; 101 PetscInt mbs = A->rmap->n, *bi = b->i, *bj = b->j; 102 MatScalar *ba = b->a; 103 PetscReal shiftnz = info->shiftamount; 104 PetscReal droptol = -1; 105 PetscBool perm_identity; 106 spbas_matrix Pattern, matrix_L, matrix_LT; 107 PetscReal mem_reduction; 108 109 PetscFunctionBegin; 110 /* Reduce memory requirements: erase values of B-matrix */ 111 PetscCall(PetscFree(ba)); 112 /* Compress (maximum) sparseness pattern of B-matrix */ 113 PetscCall(spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS, &Pattern, &mem_reduction)); 114 PetscCall(PetscFree(bi)); 115 PetscCall(PetscFree(bj)); 116 117 PetscCall(PetscInfo(NULL, " compression rate for spbas_compress_pattern %g \n", (double)mem_reduction)); 118 119 /* Make Cholesky decompositions with larger Manteuffel shifts until no more negative diagonals are found. */ 120 PetscCall(ISGetIndices(ip, &rip)); 121 PetscCall(ISGetIndices(iip, &riip)); 122 123 if (info->usedt) droptol = info->dt; 124 125 for (int ierr = NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL;) { 126 PetscBool success; 127 128 ierr = (int)spbas_incomplete_cholesky(A, rip, riip, Pattern, droptol, shiftnz, &matrix_LT, &success); 129 if (!success) { 130 shiftnz *= 1.5; 131 if (shiftnz < 1e-5) shiftnz = 1e-5; 132 PetscCall(PetscInfo(NULL, "spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%g\n", (double)shiftnz)); 133 } 134 } 135 PetscCall(spbas_delete(Pattern)); 136 137 PetscCall(PetscInfo(NULL, " memory_usage for spbas_incomplete_cholesky %g bytes per row\n", (double)(spbas_memory_requirement(matrix_LT) / (PetscReal)mbs))); 138 139 PetscCall(ISRestoreIndices(ip, &rip)); 140 PetscCall(ISRestoreIndices(iip, &riip)); 141 142 /* Convert spbas_matrix to compressed row storage */ 143 PetscCall(spbas_transpose(matrix_LT, &matrix_L)); 144 PetscCall(spbas_delete(matrix_LT)); 145 PetscCall(spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj)); 146 b->i = bi; 147 b->j = bj; 148 b->a = ba; 149 PetscCall(spbas_delete(matrix_L)); 150 151 /* Set the appropriate solution functions */ 152 PetscCall(ISIdentity(ip, &perm_identity)); 153 if (perm_identity) { 154 B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 155 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 156 B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 157 B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 158 } else { 159 B->ops->solve = MatSolve_SeqSBAIJ_1_inplace; 160 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace; 161 B->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace; 162 B->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace; 163 } 164 165 C->assembled = PETSC_TRUE; 166 C->preallocated = PETSC_TRUE; 167 168 PetscCall(PetscLogFlops(C->rmap->n)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 static PetscErrorCode MatFactorGetSolverType_seqaij_bas(Mat A, MatSolverType *type) 173 { 174 PetscFunctionBegin; 175 *type = MATSOLVERBAS; 176 PetscFunctionReturn(PETSC_SUCCESS); 177 } 178 179 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat A, MatFactorType ftype, Mat *B) 180 { 181 PetscInt n = A->rmap->n; 182 183 PetscFunctionBegin; 184 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 185 PetscCall(MatSetSizes(*B, n, n, n, n)); 186 PetscCheck(ftype == MAT_FACTOR_ICC, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 187 PetscCall(MatSetType(*B, MATSEQSBAIJ)); 188 PetscCall(MatSeqSBAIJSetPreallocation(*B, 1, MAT_SKIP_ALLOCATION, NULL)); 189 190 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqAIJ_Bas; 191 (*B)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_Bas; 192 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_bas)); 193 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_LU])); 194 PetscCall(PetscStrallocpy(MATORDERINGND, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 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(PETSC_SUCCESS); 202 } 203