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