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: }