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 Consult SPQR documentation for more information about the Common parameters 294 which correspond to the options database keys below. 295 296 Level: beginner 297 298 Note: SPQR is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 299 300 .seealso: PCQR, PCFactorSetMatSolverType(), MatSolverType 301 M*/ 302 303 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A,MatFactorType ftype,Mat *F) 304 { 305 Mat B; 306 Mat_CHOLMOD *chol; 307 PetscErrorCode ierr; 308 PetscInt m=A->rmap->n,n=A->cmap->n; 309 const char *prefix; 310 311 PetscFunctionBegin; 312 /* Create the factorization matrix F */ 313 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 314 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 315 ierr = PetscStrallocpy("spqr",&((PetscObject)B)->type_name);CHKERRQ(ierr); 316 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 317 ierr = MatSetOptionsPrefix(B,prefix);CHKERRQ(ierr); 318 ierr = MatSetUp(B);CHKERRQ(ierr); 319 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 320 321 chol->Wrap = MatWrapCholmod_SPQR_seqaij; 322 B->data = chol; 323 324 B->ops->getinfo = MatGetInfo_CHOLMOD; 325 B->ops->view = MatView_CHOLMOD; 326 B->ops->destroy = MatDestroy_CHOLMOD; 327 328 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_SPQR);CHKERRQ(ierr); 329 ierr = PetscObjectComposeFunction((PetscObject)B,"MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR);CHKERRQ(ierr); 330 331 B->factortype = MAT_FACTOR_QR; 332 B->assembled = PETSC_TRUE; 333 B->preallocated = PETSC_TRUE; 334 335 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 336 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 337 B->canuseordering = PETSC_FALSE; 338 ierr = CholmodStart(B);CHKERRQ(ierr); 339 chol->common->itype = CHOLMOD_LONG; 340 *F = B; 341 PetscFunctionReturn(0); 342 } 343 344