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