1 2 /* 3 Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1 4 5 When built with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the 6 integer type in UMFPACK, otherwise it will use int. This means 7 all integers in this file as simply declared as PetscInt. Also it means 8 that one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only] 9 10 */ 11 12 #include <../src/mat/impls/sbaij/seq/sbaij.h> 13 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 14 15 /* 16 This is a terrible hack, but it allows the error handler to retain a context. 17 Note that this hack really cannot be made both reentrant and concurrent. 18 */ 19 static Mat static_F; 20 21 static void CholmodErrorHandler(int status,const char *file,int line,const char *message) 22 { 23 PetscFunctionBegin; 24 if (status > CHOLMOD_OK) { 25 PetscCallVoid(PetscInfo(static_F,"CHOLMOD warning %d at %s:%d: %s\n",status,file,line,message)); 26 } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */ 27 PetscCallVoid(PetscInfo(static_F,"CHOLMOD OK at %s:%d: %s\n",file,line,message)); 28 } else { 29 PetscCallVoid(PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message)); 30 } 31 PetscFunctionReturnVoid(); 32 } 33 34 PetscErrorCode CholmodStart(Mat F) 35 { 36 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 37 cholmod_common *c; 38 PetscBool flg; 39 40 PetscFunctionBegin; 41 if (chol->common) PetscFunctionReturn(0); 42 PetscCall(PetscMalloc1(1,&chol->common)); 43 PetscCall(!cholmod_X_start(chol->common)); 44 45 c = chol->common; 46 c->error_handler = CholmodErrorHandler; 47 48 #define CHOLMOD_OPTION_DOUBLE(name,help) do { \ 49 PetscReal tmp = (PetscReal)c->name; \ 50 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 51 c->name = (double)tmp; \ 52 } while (0) 53 54 #define CHOLMOD_OPTION_INT(name,help) do { \ 55 PetscInt tmp = (PetscInt)c->name; \ 56 PetscCall(PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 57 c->name = (int)tmp; \ 58 } while (0) 59 60 #define CHOLMOD_OPTION_SIZE_T(name,help) do { \ 61 PetscReal tmp = (PetscInt)c->name; \ 62 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 63 PetscCheck(tmp >= 0,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"value must be positive"); \ 64 c->name = (size_t)tmp; \ 65 } while (0) 66 67 #define CHOLMOD_OPTION_BOOL(name,help) do { \ 68 PetscBool tmp = (PetscBool) !!c->name; \ 69 PetscCall(PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 70 c->name = (int)tmp; \ 71 } while (0) 72 73 PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat"); 74 CHOLMOD_OPTION_INT(nmethods,"Number of different ordering methods to try"); 75 76 #if defined(PETSC_USE_SUITESPARSE_GPU) 77 c->useGPU = 1; 78 CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0"); 79 CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes,"Maximum memory to allocate on the GPU"); 80 CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction,"Fraction of available GPU memory to allocate"); 81 #endif 82 83 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 84 chol->pack = (PetscBool)c->final_pack; 85 PetscCall(PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL)); 86 c->final_pack = (int)chol->pack; 87 88 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 89 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 90 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 91 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 92 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 93 { 94 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 95 PetscCall(PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL)); 96 } 97 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 98 CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\""); 99 CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 100 if (!c->final_asis) { 101 CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial"); 102 CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'"); 103 CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done"); 104 CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 105 } 106 { 107 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 108 PetscInt n = 3; 109 PetscCall(PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 110 PetscCheckFalse(flg && n != 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 111 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 112 } 113 { 114 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 115 PetscCall(PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 116 PetscCheckFalse(flg && n != 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 117 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 118 } 119 CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 120 CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 121 CHOLMOD_OPTION_INT(print,"Verbosity level"); 122 PetscOptionsEnd(); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc,PetscBool *valloc) 127 { 128 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)A->data; 129 PetscBool vallocin = PETSC_FALSE; 130 131 PetscFunctionBegin; 132 PetscCall(PetscMemzero(C,sizeof(*C))); 133 /* CHOLMOD uses column alignment, SBAIJ stores the upper factor, so we pass it on as a lower factor, swapping the meaning of row and column */ 134 C->nrow = (size_t)A->cmap->n; 135 C->ncol = (size_t)A->rmap->n; 136 C->nzmax = (size_t)sbaij->maxnz; 137 C->p = sbaij->i; 138 C->i = sbaij->j; 139 if (values) { 140 #if defined(PETSC_USE_COMPLEX) 141 /* we need to pass CHOLMOD the conjugate matrix */ 142 PetscScalar *v; 143 PetscInt i; 144 145 PetscCall(PetscMalloc1(sbaij->maxnz,&v)); 146 for (i = 0; i < sbaij->maxnz; i++) v[i] = PetscConj(sbaij->a[i]); 147 C->x = v; 148 vallocin = PETSC_TRUE; 149 #else 150 C->x = sbaij->a; 151 #endif 152 } 153 C->stype = -1; 154 C->itype = CHOLMOD_INT_TYPE; 155 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 156 C->dtype = CHOLMOD_DOUBLE; 157 C->sorted = 1; 158 C->packed = 1; 159 *aijalloc = PETSC_FALSE; 160 *valloc = vallocin; 161 PetscFunctionReturn(0); 162 } 163 164 #define GET_ARRAY_READ 0 165 #define GET_ARRAY_WRITE 1 166 167 PetscErrorCode VecWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 168 { 169 PetscScalar *x; 170 PetscInt n; 171 172 PetscFunctionBegin; 173 PetscCall(PetscMemzero(Y,sizeof(*Y))); 174 switch (rw) { 175 case GET_ARRAY_READ: 176 PetscCall(VecGetArrayRead(X,(const PetscScalar**)&x)); 177 break; 178 case GET_ARRAY_WRITE: 179 PetscCall(VecGetArrayWrite(X,&x)); 180 break; 181 default: 182 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 183 break; 184 } 185 PetscCall(VecGetSize(X,&n)); 186 187 Y->x = x; 188 Y->nrow = n; 189 Y->ncol = 1; 190 Y->nzmax = n; 191 Y->d = n; 192 Y->xtype = CHOLMOD_SCALAR_TYPE; 193 Y->dtype = CHOLMOD_DOUBLE; 194 PetscFunctionReturn(0); 195 } 196 197 PetscErrorCode VecUnWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 198 { 199 PetscFunctionBegin; 200 switch (rw) { 201 case GET_ARRAY_READ: 202 PetscCall(VecRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 203 break; 204 case GET_ARRAY_WRITE: 205 PetscCall(VecRestoreArrayWrite(X,(PetscScalar**)&Y->x)); 206 break; 207 default: 208 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 209 break; 210 } 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode MatDenseWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 215 { 216 PetscScalar *x; 217 PetscInt m,n,lda; 218 219 PetscFunctionBegin; 220 PetscCall(PetscMemzero(Y,sizeof(*Y))); 221 switch (rw) { 222 case GET_ARRAY_READ: 223 PetscCall(MatDenseGetArrayRead(X,(const PetscScalar**)&x)); 224 break; 225 case GET_ARRAY_WRITE: 226 PetscCall(MatDenseGetArrayWrite(X,&x)); 227 break; 228 default: 229 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 230 break; 231 } 232 PetscCall(MatDenseGetLDA(X,&lda)); 233 PetscCall(MatGetLocalSize(X,&m,&n)); 234 235 Y->x = x; 236 Y->nrow = m; 237 Y->ncol = n; 238 Y->nzmax = lda*n; 239 Y->d = lda; 240 Y->xtype = CHOLMOD_SCALAR_TYPE; 241 Y->dtype = CHOLMOD_DOUBLE; 242 PetscFunctionReturn(0); 243 } 244 245 PetscErrorCode MatDenseUnWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 246 { 247 PetscFunctionBegin; 248 switch (rw) { 249 case GET_ARRAY_READ: 250 PetscCall(MatDenseRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 251 break; 252 case GET_ARRAY_WRITE: 253 /* we don't have MatDenseRestoreArrayWrite */ 254 PetscCall(MatDenseRestoreArray(X,(PetscScalar**)&Y->x)); 255 break; 256 default: 257 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 258 break; 259 } 260 PetscFunctionReturn(0); 261 } 262 263 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 264 { 265 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 266 267 PetscFunctionBegin; 268 if (chol->spqrfact) { 269 PetscCall(!SuiteSparseQR_C_free(&chol->spqrfact, chol->common)); 270 } 271 if (chol->factor) { 272 PetscCall(!cholmod_X_free_factor(&chol->factor,chol->common)); 273 } 274 if (chol->common->itype == CHOLMOD_INT) { 275 PetscCall(!cholmod_finish(chol->common)); 276 } else { 277 PetscCall(!cholmod_l_finish(chol->common)); 278 } 279 PetscCall(PetscFree(chol->common)); 280 PetscCall(PetscFree(chol->matrix)); 281 PetscCall(PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL)); 282 PetscCall(PetscFree(F->data)); 283 PetscFunctionReturn(0); 284 } 285 286 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 287 static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat); 288 289 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 290 291 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 292 { 293 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 294 const cholmod_common *c = chol->common; 295 PetscInt i; 296 297 PetscFunctionBegin; 298 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 299 PetscCall(PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n")); 300 PetscCall(PetscViewerASCIIPushTab(viewer)); 301 PetscCall(PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE")); 302 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound)); 303 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0)); 304 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1)); 305 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2)); 306 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank)); 307 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch)); 308 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal)); 309 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis)); 310 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super)); 311 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll)); 312 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack)); 313 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic)); 314 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol)); 315 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2])); 316 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2])); 317 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper)); 318 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print)); 319 for (i=0; i<c->nmethods; i++) { 320 PetscCall(PetscViewerASCIIPrintf(viewer,"Ordering method %" PetscInt_FMT "%s:\n",i,i==c->selected ? " [SELECTED]" : "")); 321 PetscCall(PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 322 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2)); 323 } 324 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder)); 325 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis)); 326 /* Statistics */ 327 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl)); 328 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz)); 329 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz)); 330 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl)); 331 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count)); 332 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage)); 333 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse)); 334 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col)); 335 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor)); 336 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit)); 337 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl)); 338 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl)); 339 #if defined(PETSC_USE_SUITESPARSE_GPU) 340 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU)); 341 #endif 342 PetscCall(PetscViewerASCIIPopTab(viewer)); 343 PetscFunctionReturn(0); 344 } 345 346 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 347 { 348 PetscBool iascii; 349 PetscViewerFormat format; 350 351 PetscFunctionBegin; 352 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 353 if (iascii) { 354 PetscCall(PetscViewerGetFormat(viewer,&format)); 355 if (format == PETSC_VIEWER_ASCII_INFO) { 356 PetscCall(MatView_Info_CHOLMOD(F,viewer)); 357 } 358 } 359 PetscFunctionReturn(0); 360 } 361 362 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 363 { 364 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 365 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 366 367 PetscFunctionBegin; 368 static_F = F; 369 PetscCall(VecWrapCholmod(B,GET_ARRAY_READ,&cholB)); 370 PetscCall(VecWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 371 X_handle = &cholX; 372 PetscCall(!cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common)); 373 PetscCall(!cholmod_X_free_dense(&Y_handle,chol->common)); 374 PetscCall(!cholmod_X_free_dense(&E_handle,chol->common)); 375 PetscCall(VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB)); 376 PetscCall(VecUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 377 PetscCall(PetscLogFlops(4.0*chol->common->lnz)); 378 PetscFunctionReturn(0); 379 } 380 381 static PetscErrorCode MatMatSolve_CHOLMOD(Mat F,Mat B,Mat X) 382 { 383 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 384 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 385 386 PetscFunctionBegin; 387 static_F = F; 388 PetscCall(MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB)); 389 PetscCall(MatDenseWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 390 X_handle = &cholX; 391 PetscCall(!cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common)); 392 PetscCall(!cholmod_X_free_dense(&Y_handle,chol->common)); 393 PetscCall(!cholmod_X_free_dense(&E_handle,chol->common)); 394 PetscCall(MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB)); 395 PetscCall(MatDenseUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 396 PetscCall(PetscLogFlops(4.0*B->cmap->n*chol->common->lnz)); 397 PetscFunctionReturn(0); 398 } 399 400 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 401 { 402 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 403 cholmod_sparse cholA; 404 PetscBool aijalloc,valloc; 405 int err; 406 407 PetscFunctionBegin; 408 PetscCall((*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc)); 409 static_F = F; 410 err = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 411 PetscCheck(!err,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 412 PetscCheck(chol->common->status != CHOLMOD_NOT_POSDEF,PetscObjectComm((PetscObject)F),PETSC_ERR_MAT_CH_ZRPVT,"CHOLMOD detected that the matrix is not positive definite, failure at column %u",(unsigned)chol->factor->minor); 413 414 PetscCall(PetscLogFlops(chol->common->fl)); 415 if (aijalloc) PetscCall(PetscFree2(cholA.p,cholA.i)); 416 if (valloc) PetscCall(PetscFree(cholA.x)); 417 #if defined(PETSC_USE_SUITESPARSE_GPU) 418 PetscCall(PetscLogGpuTimeAdd(chol->common->CHOLMOD_GPU_GEMM_TIME + chol->common->CHOLMOD_GPU_SYRK_TIME + chol->common->CHOLMOD_GPU_TRSM_TIME + chol->common->CHOLMOD_GPU_POTRF_TIME)); 419 #endif 420 421 F->ops->solve = MatSolve_CHOLMOD; 422 F->ops->solvetranspose = MatSolve_CHOLMOD; 423 F->ops->matsolve = MatMatSolve_CHOLMOD; 424 F->ops->matsolvetranspose = MatMatSolve_CHOLMOD; 425 PetscFunctionReturn(0); 426 } 427 428 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 429 { 430 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 431 int err; 432 cholmod_sparse cholA; 433 PetscBool aijalloc,valloc; 434 PetscInt *fset = 0; 435 size_t fsize = 0; 436 437 PetscFunctionBegin; 438 PetscCall((*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc,&valloc)); 439 static_F = F; 440 if (chol->factor) { 441 err = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 442 PetscCheck(!err,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 443 } else if (perm) { 444 const PetscInt *ip; 445 PetscCall(ISGetIndices(perm,&ip)); 446 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 447 PetscCheck(chol->factor,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed using PETSc ordering with status %d",chol->common->status); 448 PetscCall(ISRestoreIndices(perm,&ip)); 449 } else { 450 chol->factor = cholmod_X_analyze(&cholA,chol->common); 451 PetscCheck(chol->factor,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed using internal ordering with status %d",chol->common->status); 452 } 453 454 if (aijalloc) PetscCall(PetscFree2(cholA.p,cholA.i)); 455 if (valloc) PetscCall(PetscFree(cholA.x)); 456 457 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 458 PetscFunctionReturn(0); 459 } 460 461 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 462 { 463 PetscFunctionBegin; 464 *type = MATSOLVERCHOLMOD; 465 PetscFunctionReturn(0); 466 } 467 468 PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F,MatInfoType flag,MatInfo *info) 469 { 470 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 471 472 PetscFunctionBegin; 473 info->block_size = 1.0; 474 info->nz_allocated = chol->common->lnz; 475 info->nz_used = chol->common->lnz; 476 info->nz_unneeded = 0.0; 477 info->assemblies = 0.0; 478 info->mallocs = 0.0; 479 info->memory = chol->common->memory_inuse; 480 info->fill_ratio_given = 0; 481 info->fill_ratio_needed = 0; 482 info->factor_mallocs = chol->common->malloc_count; 483 PetscFunctionReturn(0); 484 } 485 486 /*MC 487 MATSOLVERCHOLMOD 488 489 A matrix type providing direct solvers (Cholesky) for sequential matrices 490 via the external package CHOLMOD. 491 492 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 493 494 Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver 495 496 Consult CHOLMOD documentation for more information about the Common parameters 497 which correspond to the options database keys below. 498 499 Options Database Keys: 500 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 501 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 502 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 503 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 504 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 505 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 506 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 507 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 508 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 509 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 510 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 511 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 512 . -mat_cholmod_print <3> - Verbosity level (None) 513 - -mat_ordering_type internal - Use the ordering provided by Cholmod 514 515 Level: beginner 516 517 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 518 519 .seealso: PCCHOLESKY, PCFactorSetMatSolverType(), MatSolverType 520 M*/ 521 522 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 523 { 524 Mat B; 525 Mat_CHOLMOD *chol; 526 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 527 const char *prefix; 528 529 PetscFunctionBegin; 530 PetscCall(MatGetBlockSize(A,&bs)); 531 PetscCheck(bs == 1,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %" PetscInt_FMT,bs); 532 #if defined(PETSC_USE_COMPLEX) 533 PetscCheck(A->hermitian,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for Hermitian matrices"); 534 #endif 535 /* Create the factorization matrix F */ 536 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 537 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 538 PetscCall(PetscStrallocpy("cholmod",&((PetscObject)B)->type_name)); 539 PetscCall(MatGetOptionsPrefix(A,&prefix)); 540 PetscCall(MatSetOptionsPrefix(B,prefix)); 541 PetscCall(MatSetUp(B)); 542 PetscCall(PetscNewLog(B,&chol)); 543 544 chol->Wrap = MatWrapCholmod_seqsbaij; 545 B->data = chol; 546 547 B->ops->getinfo = MatGetInfo_CHOLMOD; 548 B->ops->view = MatView_CHOLMOD; 549 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 550 B->ops->destroy = MatDestroy_CHOLMOD; 551 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod)); 552 B->factortype = MAT_FACTOR_CHOLESKY; 553 B->assembled = PETSC_TRUE; 554 B->preallocated = PETSC_TRUE; 555 556 PetscCall(CholmodStart(B)); 557 558 PetscCall(PetscFree(B->solvertype)); 559 PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype)); 560 B->canuseordering = PETSC_TRUE; 561 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_CHOLESKY])); 562 *F = B; 563 PetscFunctionReturn(0); 564 } 565