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 = PetscInfo4(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 = PetscInfo3(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 if (tmp < 0) SETERRQ(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 if (flg && n != 3) SETERRQ(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 if (flg && n != 3) SETERRQ(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 static 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 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D 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 static 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 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D not handled",rw); 216 break; 217 } 218 PetscFunctionReturn(0); 219 } 220 221 static 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 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D 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 static 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 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D 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 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 280 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 281 ierr = PetscFree(chol->common);CHKERRQ(ierr); 282 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 283 ierr = PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr); 284 ierr = PetscFree(F->data);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 289 static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat); 290 291 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 292 293 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 294 { 295 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 296 const cholmod_common *c = chol->common; 297 PetscErrorCode ierr; 298 PetscInt i; 299 300 PetscFunctionBegin; 301 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 302 ierr = PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");CHKERRQ(ierr); 303 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 304 ierr = PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");CHKERRQ(ierr); 305 ierr = PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound);CHKERRQ(ierr); 306 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0);CHKERRQ(ierr); 307 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1);CHKERRQ(ierr); 308 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2);CHKERRQ(ierr); 309 ierr = PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank);CHKERRQ(ierr); 310 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);CHKERRQ(ierr); 311 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal);CHKERRQ(ierr); 312 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis);CHKERRQ(ierr); 313 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super);CHKERRQ(ierr); 314 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll);CHKERRQ(ierr); 315 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack);CHKERRQ(ierr); 316 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic);CHKERRQ(ierr); 317 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol);CHKERRQ(ierr); 318 ierr = PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print);CHKERRQ(ierr); 322 for (i=0; i<c->nmethods; i++) { 323 ierr = PetscViewerASCIIPrintf(viewer,"Ordering method %D%s:\n",i,i==c->selected ? " [SELECTED]" : "");CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 325 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2);CHKERRQ(ierr); 326 } 327 ierr = PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);CHKERRQ(ierr); 329 /* Statistics */ 330 ierr = PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz);CHKERRQ(ierr); 332 ierr = PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz);CHKERRQ(ierr); 333 ierr = PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl);CHKERRQ(ierr); 334 ierr = PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count);CHKERRQ(ierr); 335 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);CHKERRQ(ierr); 341 ierr = PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);CHKERRQ(ierr); 342 #if defined(PETSC_USE_SUITESPARSE_GPU) 343 ierr = PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);CHKERRQ(ierr); 344 #endif 345 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 350 { 351 PetscErrorCode ierr; 352 PetscBool iascii; 353 PetscViewerFormat format; 354 355 PetscFunctionBegin; 356 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 357 if (iascii) { 358 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 359 if (format == PETSC_VIEWER_ASCII_INFO) { 360 ierr = MatView_Info_CHOLMOD(F,viewer);CHKERRQ(ierr); 361 } 362 } 363 PetscFunctionReturn(0); 364 } 365 366 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 367 { 368 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 369 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 static_F = F; 374 ierr = VecWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 375 ierr = VecWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 376 X_handle = &cholX; 377 ierr = !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);CHKERRQ(ierr); 378 ierr = !cholmod_X_free_dense(&Y_handle,chol->common);CHKERRQ(ierr); 379 ierr = !cholmod_X_free_dense(&E_handle,chol->common);CHKERRQ(ierr); 380 ierr = VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 381 ierr = VecUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 382 PetscFunctionReturn(0); 383 } 384 385 static PetscErrorCode MatMatSolve_CHOLMOD(Mat F,Mat B,Mat X) 386 { 387 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 388 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 389 PetscErrorCode ierr; 390 391 PetscFunctionBegin; 392 static_F = F; 393 ierr = MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 394 ierr = MatDenseWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 395 X_handle = &cholX; 396 ierr = !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);CHKERRQ(ierr); 397 ierr = !cholmod_X_free_dense(&Y_handle,chol->common);CHKERRQ(ierr); 398 ierr = !cholmod_X_free_dense(&E_handle,chol->common);CHKERRQ(ierr); 399 ierr = MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 400 ierr = MatDenseUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 405 { 406 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 407 cholmod_sparse cholA; 408 PetscBool aijalloc,valloc; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 413 static_F = F; 414 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 415 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 416 if (chol->common->status == CHOLMOD_NOT_POSDEF) SETERRQ1(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); 417 418 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 419 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 420 #if defined(PETSC_USE_SUITESPARSE_GPU) 421 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); 422 #endif 423 424 F->ops->solve = MatSolve_CHOLMOD; 425 F->ops->solvetranspose = MatSolve_CHOLMOD; 426 F->ops->matsolve = MatMatSolve_CHOLMOD; 427 F->ops->matsolvetranspose = MatMatSolve_CHOLMOD; 428 PetscFunctionReturn(0); 429 } 430 431 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 432 { 433 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 434 PetscErrorCode ierr; 435 cholmod_sparse cholA; 436 PetscBool aijalloc,valloc; 437 PetscInt *fset = 0; 438 size_t fsize = 0; 439 440 PetscFunctionBegin; 441 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 442 static_F = F; 443 if (chol->factor) { 444 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 445 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 446 } else if (perm) { 447 const PetscInt *ip; 448 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 449 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 450 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 451 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 452 } else { 453 chol->factor = cholmod_X_analyze(&cholA,chol->common); 454 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 455 } 456 457 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 458 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 459 460 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 461 PetscFunctionReturn(0); 462 } 463 464 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 465 { 466 PetscFunctionBegin; 467 *type = MATSOLVERCHOLMOD; 468 PetscFunctionReturn(0); 469 } 470 471 PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F,MatInfoType flag,MatInfo *info) 472 { 473 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 474 475 PetscFunctionBegin; 476 info->block_size = 1.0; 477 info->nz_allocated = chol->common->lnz; 478 info->nz_used = chol->common->lnz; 479 info->nz_unneeded = 0.0; 480 info->assemblies = 0.0; 481 info->mallocs = 0.0; 482 info->memory = chol->common->memory_inuse; 483 info->fill_ratio_given = 0; 484 info->fill_ratio_needed = 0; 485 info->factor_mallocs = chol->common->malloc_count; 486 PetscFunctionReturn(0); 487 } 488 489 /*MC 490 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 491 via the external package CHOLMOD. 492 493 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 494 495 Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver 496 497 Consult CHOLMOD documentation for more information about the Common parameters 498 which correspond to the options database keys below. 499 500 Options Database Keys: 501 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 502 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 503 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 504 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 505 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 506 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 507 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 508 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 509 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 510 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 511 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 512 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 513 - -mat_cholmod_print <3> - Verbosity level (None) 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 PetscErrorCode ierr; 527 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 528 const char *prefix; 529 530 PetscFunctionBegin; 531 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 532 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 533 #if defined(PETSC_USE_COMPLEX) 534 if (!A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for hermitian matrices"); 535 #endif 536 /* Create the factorization matrix F */ 537 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 538 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 539 ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr); 540 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 541 ierr = MatSetOptionsPrefix(B,prefix);CHKERRQ(ierr); 542 ierr = MatSetUp(B);CHKERRQ(ierr); 543 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 544 545 chol->Wrap = MatWrapCholmod_seqsbaij; 546 B->data = chol; 547 548 B->ops->getinfo = MatGetInfo_CHOLMOD; 549 B->ops->view = MatView_CHOLMOD; 550 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 551 B->ops->destroy = MatDestroy_CHOLMOD; 552 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);CHKERRQ(ierr); 553 B->factortype = MAT_FACTOR_CHOLESKY; 554 B->assembled = PETSC_TRUE; 555 B->preallocated = PETSC_TRUE; 556 557 ierr = CholmodStart(B);CHKERRQ(ierr); 558 559 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 560 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 561 562 *F = B; 563 PetscFunctionReturn(0); 564 } 565