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