1 2 #include <petscsys.h> 3 #include <../src/mat/impls/aij/seq/aij.h> 4 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 5 6 EXTERN_C_BEGIN 7 #include <SuiteSparseQR_C.h> 8 EXTERN_C_END 9 10 static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc,PetscBool *valloc) 11 { 12 Mat_SeqAIJ *aij; 13 Mat AT; 14 const PetscScalar *aa; 15 PetscScalar *ca; 16 const PetscInt *ai, *aj; 17 PetscInt n = A->cmap->n, i,j,k,nz; 18 SuiteSparse_long *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */ 19 PetscBool vain = PETSC_FALSE,flg; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg);CHKERRQ(ierr); 24 if (flg) { 25 ierr = MatNormalGetMat(A, &A);CHKERRQ(ierr); 26 } 27 /* cholmod_sparse is compressed sparse column */ 28 ierr = MatGetOption(A, MAT_SYMMETRIC, &flg);CHKERRQ(ierr); 29 if (flg) { 30 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 31 AT = A; 32 } else { 33 ierr = MatTranspose(A, MAT_INITIAL_MATRIX, &AT);CHKERRQ(ierr); 34 } 35 aij = (Mat_SeqAIJ*)AT->data; 36 ai = aij->j; 37 aj = aij->i; 38 for (j=0,nz=0; j<n; j++) nz += aj[j+1] - aj[j]; 39 ierr = PetscMalloc2(n+1,&cj,nz,&ci);CHKERRQ(ierr); 40 if (values) { 41 vain = PETSC_TRUE; 42 ierr = PetscMalloc1(nz,&ca);CHKERRQ(ierr); 43 ierr = MatSeqAIJGetArrayRead(AT,&aa);CHKERRQ(ierr); 44 } 45 for (j=0,k=0; j<n; j++) { 46 cj[j] = k; 47 for (i=aj[j]; i<aj[j+1]; i++,k++) { 48 ci[k] = ai[i]; 49 if (values) ca[k] = aa[i]; 50 } 51 } 52 cj[j] = k; 53 *aijalloc = PETSC_TRUE; 54 *valloc = vain; 55 if (values) { 56 ierr = MatSeqAIJRestoreArrayRead(AT,&aa);CHKERRQ(ierr); 57 } 58 59 ierr = PetscMemzero(C,sizeof(*C));CHKERRQ(ierr); 60 61 C->nrow = (size_t)AT->cmap->n; 62 C->ncol = (size_t)AT->rmap->n; 63 C->nzmax = (size_t)nz; 64 C->p = cj; 65 C->i = ci; 66 C->x = values ? ca : 0; 67 C->stype = 0; 68 C->itype = CHOLMOD_LONG; 69 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 70 C->dtype = CHOLMOD_DOUBLE; 71 C->sorted = 1; 72 C->packed = 1; 73 74 ierr = MatDestroy(&AT);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A,MatSolverType *type) 79 { 80 PetscFunctionBegin; 81 *type = MATSOLVERSPQR; 82 PetscFunctionReturn(0); 83 } 84 85 #define GET_ARRAY_READ 0 86 #define GET_ARRAY_WRITE 1 87 88 static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle) 89 { 90 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 91 cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 if (!chol->normal) { 96 QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common); 97 if (!QTB_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed"); 98 Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common); 99 if (!Y_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 100 } else { 101 Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common); 102 if (!Z_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 103 Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common); 104 if (!Y_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 105 ierr = !cholmod_l_free_dense(&Z_handle, chol->common);CHKERRQ(ierr); 106 } 107 *_Y_handle = Y_handle; 108 ierr = !cholmod_l_free_dense(&QTB_handle, chol->common);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode MatSolve_SPQR(Mat F,Vec B,Vec X) 113 { 114 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 115 cholmod_dense cholB,*Y_handle = NULL; 116 PetscInt n; 117 PetscScalar *v; 118 PetscErrorCode ierr; 119 120 PetscFunctionBegin; 121 ierr = VecWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 122 ierr = MatSolve_SPQR_Internal(F, &cholB, &Y_handle);CHKERRQ(ierr); 123 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr); 124 ierr = VecGetArrayWrite(X, &v);CHKERRQ(ierr); 125 ierr = PetscArraycpy(v, (PetscScalar *) (Y_handle->x), n);CHKERRQ(ierr); 126 ierr = VecRestoreArrayWrite(X, &v);CHKERRQ(ierr); 127 ierr = !cholmod_l_free_dense(&Y_handle, chol->common);CHKERRQ(ierr); 128 ierr = VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 132 static PetscErrorCode MatMatSolve_SPQR(Mat F,Mat B,Mat X) 133 { 134 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 135 cholmod_dense cholB,*Y_handle = NULL; 136 PetscScalar *v; 137 PetscInt lda; 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 ierr = MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 142 ierr = MatSolve_SPQR_Internal(F, &cholB, &Y_handle);CHKERRQ(ierr); 143 ierr = MatDenseGetArrayWrite(X, &v);CHKERRQ(ierr); 144 ierr = MatDenseGetLDA(X, &lda);CHKERRQ(ierr); 145 if ((size_t) lda == Y_handle->d) { 146 ierr = PetscArraycpy(v, (PetscScalar *) (Y_handle->x), lda * Y_handle->ncol);CHKERRQ(ierr); 147 } else { 148 for (size_t j = 0; j < Y_handle->ncol; j++) { 149 ierr = PetscArraycpy(&v[j*lda], &(((PetscScalar *) Y_handle->x)[j*Y_handle->d]), Y_handle->nrow);CHKERRQ(ierr); 150 } 151 } 152 ierr = MatDenseRestoreArrayWrite(X, &v);CHKERRQ(ierr); 153 ierr = !cholmod_l_free_dense(&Y_handle, chol->common);CHKERRQ(ierr); 154 ierr = MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle) 159 { 160 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 161 cholmod_dense *Y_handle = NULL, *RTB_handle = NULL; 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common); 166 if (!RTB_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 167 Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common); 168 if (!Y_handle) SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed"); 169 *_Y_handle = Y_handle; 170 ierr = !cholmod_l_free_dense(&RTB_handle, chol->common);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 static PetscErrorCode MatSolveTranspose_SPQR(Mat F,Vec B,Vec X) 175 { 176 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 177 cholmod_dense cholB,*Y_handle = NULL; 178 PetscInt n; 179 PetscScalar *v; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 ierr = VecWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 184 ierr = MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);CHKERRQ(ierr); 185 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr); 186 ierr = VecGetArrayWrite(X, &v);CHKERRQ(ierr); 187 ierr = PetscArraycpy(v, (PetscScalar *) Y_handle->x, n);CHKERRQ(ierr); 188 ierr = VecRestoreArrayWrite(X, &v);CHKERRQ(ierr); 189 ierr = !cholmod_l_free_dense(&Y_handle, chol->common);CHKERRQ(ierr); 190 ierr = VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F,Mat B,Mat X) 195 { 196 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 197 cholmod_dense cholB,*Y_handle = NULL; 198 PetscScalar *v; 199 PetscInt lda; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 204 ierr = MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);CHKERRQ(ierr); 205 ierr = MatDenseGetArrayWrite(X, &v);CHKERRQ(ierr); 206 ierr = MatDenseGetLDA(X, &lda);CHKERRQ(ierr); 207 if ((size_t) lda == Y_handle->d) { 208 ierr = PetscArraycpy(v, (PetscScalar *) Y_handle->x, lda * Y_handle->ncol);CHKERRQ(ierr); 209 } else { 210 for (size_t j = 0; j < Y_handle->ncol; j++) { 211 ierr = PetscArraycpy(&v[j*lda], &(((PetscScalar *) Y_handle->x)[j*Y_handle->d]), Y_handle->nrow);CHKERRQ(ierr); 212 } 213 } 214 ierr = MatDenseRestoreArrayWrite(X, &v);CHKERRQ(ierr); 215 ierr = !cholmod_l_free_dense(&Y_handle, chol->common);CHKERRQ(ierr); 216 ierr = MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F,Mat A,const MatFactorInfo *info) 221 { 222 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 223 cholmod_sparse cholA; 224 PetscBool aijalloc,valloc,flg; 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 ierr = PetscObjectTypeCompare((PetscObject)A,MATNORMAL,&flg);CHKERRQ(ierr); 229 chol->normal = flg; 230 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 231 ierr = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common); 232 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"SPQR factorization failed with status %d",chol->common->status); 233 234 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 235 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 236 237 F->ops->solve = MatSolve_SPQR; 238 F->ops->matsolve = MatMatSolve_SPQR; 239 if (chol->normal) { 240 F->ops->solvetranspose = MatSolve_SPQR; 241 F->ops->matsolvetranspose = MatMatSolve_SPQR; 242 } else if (A->cmap->n == A->rmap->n) { 243 F->ops->solvetranspose = MatSolveTranspose_SPQR; 244 F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR; 245 } 246 PetscFunctionReturn(0); 247 } 248 249 PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F,Mat A,IS perm,const MatFactorInfo *info) 250 { 251 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 252 PetscErrorCode ierr; 253 cholmod_sparse cholA; 254 PetscBool aijalloc,valloc,flg; 255 256 PetscFunctionBegin; 257 ierr = PetscObjectTypeCompare((PetscObject)A,MATNORMAL,&flg);CHKERRQ(ierr); 258 chol->normal = flg; 259 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 260 if (PetscDefined(USE_DEBUG)) { 261 ierr = !cholmod_l_check_sparse(&cholA, chol->common);CHKERRQ(ierr); 262 } 263 if (chol->spqrfact) { 264 ierr = !SuiteSparseQR_C_free(&chol->spqrfact, chol->common);CHKERRQ(ierr); 265 } 266 chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common); 267 if (!chol->spqrfact) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed using internal ordering with status %d",chol->common->status); 268 269 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 270 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 271 272 ierr = PetscObjectComposeFunction((PetscObject)F,"MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR);CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 /*MC 277 MATSOLVERSPQR 278 279 A matrix type providing direct solvers (QR factorizations) for sequential matrices 280 via the external package SPQR. 281 282 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 283 284 TODO: Use -pc_type qr -pc_factor_mat_solver_type spqr to use this direct solver (TODO: PCQR) 285 286 Consult SPQR documentation for more information about the Common parameters 287 which correspond to the options database keys below. 288 289 Level: beginner 290 291 Note: SPQR is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 292 293 .seealso: PCQR, PCFactorSetMatSolverType(), MatSolverType 294 M*/ 295 296 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A,MatFactorType ftype,Mat *F) 297 { 298 Mat B; 299 Mat_CHOLMOD *chol; 300 PetscErrorCode ierr; 301 PetscInt m=A->rmap->n,n=A->cmap->n; 302 const char *prefix; 303 304 PetscFunctionBegin; 305 /* Create the factorization matrix F */ 306 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 307 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 308 ierr = PetscStrallocpy("spqr",&((PetscObject)B)->type_name);CHKERRQ(ierr); 309 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 310 ierr = MatSetOptionsPrefix(B,prefix);CHKERRQ(ierr); 311 ierr = MatSetUp(B);CHKERRQ(ierr); 312 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 313 314 chol->Wrap = MatWrapCholmod_SPQR_seqaij; 315 B->data = chol; 316 317 B->ops->getinfo = MatGetInfo_CHOLMOD; 318 B->ops->view = MatView_CHOLMOD; 319 B->ops->destroy = MatDestroy_CHOLMOD; 320 321 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_SPQR);CHKERRQ(ierr); 322 ierr = PetscObjectComposeFunction((PetscObject)B,"MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR);CHKERRQ(ierr); 323 324 B->factortype = MAT_FACTOR_QR; 325 B->assembled = PETSC_TRUE; 326 B->preallocated = PETSC_TRUE; 327 328 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 329 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 330 B->canuseordering = PETSC_FALSE; 331 ierr = CholmodStart(B);CHKERRQ(ierr); 332 chol->common->itype = CHOLMOD_LONG; 333 *F = B; 334 PetscFunctionReturn(0); 335 } 336 337