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