Actual source code: aijspqr.c


  2: #include <petscsys.h>
  3: #include <../src/mat/impls/aij/seq/aij.h>
  4: #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>

  6: EXTERN_C_BEGIN
  7: #include <SuiteSparseQR_C.h>
  8: EXTERN_C_END

 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;

 23:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg);
 24:   if (flg) {
 25:     MatNormalHermitianGetMat(A, &A);
 26:   } else if (!PetscDefined(USE_COMPLEX)) {
 27:     PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg);
 28:     if (flg) {
 29:       MatNormalGetMat(A, &A);
 30:     }
 31:   }
 32:   /* cholmod_sparse is compressed sparse column */
 33:   MatGetOption(A, MAT_SYMMETRIC, &flg);
 34:   if (flg) {
 35:     PetscObjectReference((PetscObject)A);
 36:     AT = A;
 37:   } else {
 38:     MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
 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:   PetscMalloc2(n+1,&cj,nz,&ci);
 45:   if (values) {
 46:     vain = PETSC_TRUE;
 47:     PetscMalloc1(nz,&ca);
 48:     MatSeqAIJGetArrayRead(AT,&aa);
 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:     MatSeqAIJRestoreArrayRead(AT,&aa);
 62:   }

 64:   PetscMemzero(C,sizeof(*C));

 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;

 79:   MatDestroy(&AT);
 80:   return(0);
 81: }

 83: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A,MatSolverType *type)
 84: {
 86:   *type = MATSOLVERSPQR;
 87:   return(0);
 88: }

 90: #define GET_ARRAY_READ 0
 91: #define GET_ARRAY_WRITE 1

 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;

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:     !cholmod_l_free_dense(&Z_handle, chol->common);
111:   }
112:   *_Y_handle = Y_handle;
113:   !cholmod_l_free_dense(&QTB_handle, chol->common);
114:   return(0);
115: }

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;

126:   VecWrapCholmod(B,GET_ARRAY_READ,&cholB);
127:   MatSolve_SPQR_Internal(F, &cholB, &Y_handle);
128:   VecGetLocalSize(X, &n);
129:   VecGetArrayWrite(X, &v);
130:   PetscArraycpy(v, (PetscScalar *) (Y_handle->x), n);
131:   VecRestoreArrayWrite(X, &v);
132:   !cholmod_l_free_dense(&Y_handle, chol->common);
133:   VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
134:   return(0);
135: }

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;

146:   MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);
147:   MatSolve_SPQR_Internal(F, &cholB, &Y_handle);
148:   MatDenseGetArrayWrite(X, &v);
149:   MatDenseGetLDA(X, &lda);
150:   if ((size_t) lda == Y_handle->d) {
151:     PetscArraycpy(v, (PetscScalar *) (Y_handle->x), lda * Y_handle->ncol);
152:   } else {
153:     for (size_t j = 0; j < Y_handle->ncol; j++) {
154:       PetscArraycpy(&v[j*lda], &(((PetscScalar *) Y_handle->x)[j*Y_handle->d]), Y_handle->nrow);
155:     }
156:   }
157:   MatDenseRestoreArrayWrite(X, &v);
158:   !cholmod_l_free_dense(&Y_handle, chol->common);
159:   MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
160:   return(0);
161: }

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;

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:   !cholmod_l_free_dense(&RTB_handle, chol->common);
176:   return(0);
177: }

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;

188:   VecWrapCholmod(B,GET_ARRAY_READ,&cholB);
189:   MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);
190:   VecGetLocalSize(X, &n);
191:   VecGetArrayWrite(X, &v);
192:   PetscArraycpy(v, (PetscScalar *) Y_handle->x, n);
193:   VecRestoreArrayWrite(X, &v);
194:   !cholmod_l_free_dense(&Y_handle, chol->common);
195:   VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
196:   return(0);
197: }

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;

208:   MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);
209:   MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);
210:   MatDenseGetArrayWrite(X, &v);
211:   MatDenseGetLDA(X, &lda);
212:   if ((size_t) lda == Y_handle->d) {
213:     PetscArraycpy(v, (PetscScalar *) Y_handle->x, lda * Y_handle->ncol);
214:   } else {
215:     for (size_t j = 0; j < Y_handle->ncol; j++) {
216:       PetscArraycpy(&v[j*lda], &(((PetscScalar *) Y_handle->x)[j*Y_handle->d]), Y_handle->nrow);
217:     }
218:   }
219:   MatDenseRestoreArrayWrite(X, &v);
220:   !cholmod_l_free_dense(&Y_handle, chol->common);
221:   MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
222:   return(0);
223: }

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;

233:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal);
234:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) {
235:     PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal);
236:   }
237:   (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);
238:   !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);

241:   if (aijalloc) {PetscFree2(cholA.p,cholA.i);}
242:   if (valloc) {PetscFree(cholA.x);}

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:   return(0);
254: }

256: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F,Mat A,IS perm,const MatFactorInfo *info)
257: {
258:   Mat_CHOLMOD    *chol = (Mat_CHOLMOD*)F->data;
260:   cholmod_sparse cholA;
261:   PetscBool      aijalloc,valloc;

264:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal);
265:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) {
266:     PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal);
267:   }
268:   (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);
269:   if (PetscDefined(USE_DEBUG)) {
270:     !cholmod_l_check_sparse(&cholA, chol->common);
271:   }
272:   if (chol->spqrfact) {
273:     !SuiteSparseQR_C_free(&chol->spqrfact, chol->common);
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);

278:   if (aijalloc) {PetscFree2(cholA.p,cholA.i);}
279:   if (valloc) {PetscFree(cholA.x);}

281:   PetscObjectComposeFunction((PetscObject)F,"MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR);
282:   return(0);
283: }

285: /*MC
286:   MATSOLVERSPQR

288:   A matrix type providing direct solvers (QR factorizations) for sequential matrices
289:   via the external package SPQR.

291:   Use ./configure --download-suitesparse to install PETSc to use CHOLMOD

293:   Consult SPQR documentation for more information about the Common parameters
294:   which correspond to the options database keys below.

296:    Level: beginner

298:    Note: SPQR is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

300: .seealso: PCQR, PCFactorSetMatSolverType(), MatSolverType
301: M*/

303: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A,MatFactorType ftype,Mat *F)
304: {
305:   Mat            B;
306:   Mat_CHOLMOD    *chol;
308:   PetscInt       m=A->rmap->n,n=A->cmap->n;
309:   const char     *prefix;

312:   /* Create the factorization matrix F */
313:   MatCreate(PetscObjectComm((PetscObject)A),&B);
314:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
315:   PetscStrallocpy("spqr",&((PetscObject)B)->type_name);
316:   MatGetOptionsPrefix(A,&prefix);
317:   MatSetOptionsPrefix(B,prefix);
318:   MatSetUp(B);
319:   PetscNewLog(B,&chol);

321:   chol->Wrap = MatWrapCholmod_SPQR_seqaij;
322:   B->data    = chol;

324:   B->ops->getinfo = MatGetInfo_CHOLMOD;
325:   B->ops->view    = MatView_CHOLMOD;
326:   B->ops->destroy = MatDestroy_CHOLMOD;

328:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_SPQR);
329:   PetscObjectComposeFunction((PetscObject)B,"MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR);

331:   B->factortype   = MAT_FACTOR_QR;
332:   B->assembled    = PETSC_TRUE;
333:   B->preallocated = PETSC_TRUE;

335:   PetscFree(B->solvertype);
336:   PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);
337:   B->canuseordering = PETSC_FALSE;
338:   CholmodStart(B);
339:   chol->common->itype = CHOLMOD_LONG;
340:   *F   = B;
341:   return(0);
342: }