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;

 21:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg);
 22:   if (flg) {
 23:     MatNormalHermitianGetMat(A, &A);
 24:   } else if (!PetscDefined(USE_COMPLEX)) {
 25:     PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg);
 26:     if (flg) MatNormalGetMat(A, &A);
 27:   }
 28:   /* cholmod_sparse is compressed sparse column */
 29:   MatIsSymmetric(A, 0.0, &flg);
 30:   if (flg) {
 31:     PetscObjectReference((PetscObject)A);
 32:     AT = A;
 33:   } else {
 34:     MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
 35:   }
 36:   aij = (Mat_SeqAIJ *)AT->data;
 37:   ai  = aij->j;
 38:   aj  = aij->i;
 39:   for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
 40:   PetscMalloc2(n + 1, &cj, nz, &ci);
 41:   if (values) {
 42:     vain = PETSC_TRUE;
 43:     PetscMalloc1(nz, &ca);
 44:     MatSeqAIJGetArrayRead(AT, &aa);
 45:   }
 46:   for (j = 0, k = 0; j < n; j++) {
 47:     cj[j] = k;
 48:     for (i = aj[j]; i < aj[j + 1]; i++, k++) {
 49:       ci[k] = ai[i];
 50:       if (values) ca[k] = aa[i];
 51:     }
 52:   }
 53:   cj[j]     = k;
 54:   *aijalloc = PETSC_TRUE;
 55:   *valloc   = vain;
 56:   if (values) MatSeqAIJRestoreArrayRead(AT, &aa);

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

 60:   C->nrow   = (size_t)AT->cmap->n;
 61:   C->ncol   = (size_t)AT->rmap->n;
 62:   C->nzmax  = (size_t)nz;
 63:   C->p      = cj;
 64:   C->i      = ci;
 65:   C->x      = values ? ca : 0;
 66:   C->stype  = 0;
 67:   C->itype  = CHOLMOD_LONG;
 68:   C->xtype  = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
 69:   C->dtype  = CHOLMOD_DOUBLE;
 70:   C->sorted = 1;
 71:   C->packed = 1;

 73:   MatDestroy(&AT);
 74:   return 0;
 75: }

 77: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
 78: {
 79:   *type = MATSOLVERSPQR;
 80:   return 0;
 81: }

 83: #define GET_ARRAY_READ  0
 84: #define GET_ARRAY_WRITE 1

 86: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
 87: {
 88:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
 89:   cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;

 91:   if (!chol->normal) {
 92:     QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
 94:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
 96:   } else {
 97:     Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
 99:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
101:     !cholmod_l_free_dense(&Z_handle, chol->common);
102:   }
103:   *_Y_handle = Y_handle;
104:   !cholmod_l_free_dense(&QTB_handle, chol->common);
105:   return 0;
106: }

108: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
109: {
110:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
111:   cholmod_dense cholB, *Y_handle = NULL;
112:   PetscInt      n;
113:   PetscScalar  *v;

115:   VecWrapCholmod(B, GET_ARRAY_READ, &cholB);
116:   MatSolve_SPQR_Internal(F, &cholB, &Y_handle);
117:   VecGetLocalSize(X, &n);
118:   VecGetArrayWrite(X, &v);
119:   PetscArraycpy(v, (PetscScalar *)(Y_handle->x), n);
120:   VecRestoreArrayWrite(X, &v);
121:   !cholmod_l_free_dense(&Y_handle, chol->common);
122:   VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB);
123:   return 0;
124: }

126: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
127: {
128:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
129:   cholmod_dense cholB, *Y_handle = NULL;
130:   PetscScalar  *v;
131:   PetscInt      lda;

133:   MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB);
134:   MatSolve_SPQR_Internal(F, &cholB, &Y_handle);
135:   MatDenseGetArrayWrite(X, &v);
136:   MatDenseGetLDA(X, &lda);
137:   if ((size_t)lda == Y_handle->d) {
138:     PetscArraycpy(v, (PetscScalar *)(Y_handle->x), lda * Y_handle->ncol);
139:   } else {
140:     for (size_t j = 0; j < Y_handle->ncol; j++) PetscArraycpy(&v[j * lda], &(((PetscScalar *)Y_handle->x)[j * Y_handle->d]), Y_handle->nrow);
141:   }
142:   MatDenseRestoreArrayWrite(X, &v);
143:   !cholmod_l_free_dense(&Y_handle, chol->common);
144:   MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB);
145:   return 0;
146: }

148: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
149: {
150:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
151:   cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;

153:   RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
155:   Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
157:   *_Y_handle = Y_handle;
158:   !cholmod_l_free_dense(&RTB_handle, chol->common);
159:   return 0;
160: }

162: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
163: {
164:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
165:   cholmod_dense cholB, *Y_handle = NULL;
166:   PetscInt      n;
167:   PetscScalar  *v;

169:   VecWrapCholmod(B, GET_ARRAY_READ, &cholB);
170:   MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);
171:   VecGetLocalSize(X, &n);
172:   VecGetArrayWrite(X, &v);
173:   PetscArraycpy(v, (PetscScalar *)Y_handle->x, n);
174:   VecRestoreArrayWrite(X, &v);
175:   !cholmod_l_free_dense(&Y_handle, chol->common);
176:   VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB);
177:   return 0;
178: }

180: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
181: {
182:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
183:   cholmod_dense cholB, *Y_handle = NULL;
184:   PetscScalar  *v;
185:   PetscInt      lda;

187:   MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB);
188:   MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle);
189:   MatDenseGetArrayWrite(X, &v);
190:   MatDenseGetLDA(X, &lda);
191:   if ((size_t)lda == Y_handle->d) {
192:     PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol);
193:   } else {
194:     for (size_t j = 0; j < Y_handle->ncol; j++) PetscArraycpy(&v[j * lda], &(((PetscScalar *)Y_handle->x)[j * Y_handle->d]), Y_handle->nrow);
195:   }
196:   MatDenseRestoreArrayWrite(X, &v);
197:   !cholmod_l_free_dense(&Y_handle, chol->common);
198:   MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB);
199:   return 0;
200: }

202: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
203: {
204:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
205:   cholmod_sparse cholA;
206:   PetscBool      aijalloc, valloc;
207:   int            err;

209:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal);
210:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal);
211:   (*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc);
212:   err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);

215:   if (aijalloc) PetscFree2(cholA.p, cholA.i);
216:   if (valloc) PetscFree(cholA.x);

218:   F->ops->solve    = MatSolve_SPQR;
219:   F->ops->matsolve = MatMatSolve_SPQR;
220:   if (chol->normal) {
221:     F->ops->solvetranspose    = MatSolve_SPQR;
222:     F->ops->matsolvetranspose = MatMatSolve_SPQR;
223:   } else if (A->cmap->n == A->rmap->n) {
224:     F->ops->solvetranspose    = MatSolveTranspose_SPQR;
225:     F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
226:   }
227:   return 0;
228: }

230: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
231: {
232:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
233:   cholmod_sparse cholA;
234:   PetscBool      aijalloc, valloc;

236:   PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal);
237:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal);
238:   (*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc);
239:   if (PetscDefined(USE_DEBUG)) !cholmod_l_check_sparse(&cholA, chol->common);
240:   if (chol->spqrfact) !SuiteSparseQR_C_free(&chol->spqrfact, chol->common);
241:   chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);

244:   if (aijalloc) PetscFree2(cholA.p, cholA.i);
245:   if (valloc) PetscFree(cholA.x);

247:   PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR);
248:   return 0;
249: }

251: /*MC
252:   MATSOLVERSPQR

254:   A matrix type providing direct solvers, QR factorizations, for sequential matrices
255:   via the external package SPQR.

257:   Use ./configure --download-suitesparse to install PETSc to use SPQR

259:   Consult SPQR documentation for more information about the common parameters
260:   which correspond to the options database keys below.

262:    Level: beginner

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

267: .seealso: `PCQR`, `PCFactorSetMatSolverType()`, `MatSolverType`
268: M*/

270: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F)
271: {
272:   Mat          B;
273:   Mat_CHOLMOD *chol;
274:   PetscInt     m = A->rmap->n, n = A->cmap->n;
275:   const char  *prefix;

277:   /* Create the factorization matrix F */
278:   MatCreate(PetscObjectComm((PetscObject)A), &B);
279:   MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n);
280:   PetscStrallocpy("spqr", &((PetscObject)B)->type_name);
281:   MatGetOptionsPrefix(A, &prefix);
282:   MatSetOptionsPrefix(B, prefix);
283:   MatSetUp(B);
284:   PetscNew(&chol);

286:   chol->Wrap = MatWrapCholmod_SPQR_seqaij;
287:   B->data    = chol;

289:   B->ops->getinfo = MatGetInfo_CHOLMOD;
290:   B->ops->view    = MatView_CHOLMOD;
291:   B->ops->destroy = MatDestroy_CHOLMOD;

293:   PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR);
294:   PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR);

296:   B->factortype   = MAT_FACTOR_QR;
297:   B->assembled    = PETSC_TRUE;
298:   B->preallocated = PETSC_TRUE;

300:   PetscFree(B->solvertype);
301:   PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype);
302:   B->canuseordering = PETSC_FALSE;
303:   CholmodStart(B);
304:   chol->common->itype = CHOLMOD_LONG;
305:   *F                  = B;
306:   return 0;
307: }