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