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