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