Actual source code: dense.c
2: /*
3: Defines the basic matrix operations for sequential dense.
4: */
6: #include <../src/mat/impls/dense/seq/dense.h>
7: #include <petscblaslapack.h>
9: #include <../src/mat/impls/aij/seq/aij.h>
11: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
12: {
13: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
14: PetscInt j, k, n = A->rmap->n;
15: PetscScalar *v;
18: MatDenseGetArray(A,&v);
19: if (!hermitian) {
20: for (k=0;k<n;k++) {
21: for (j=k;j<n;j++) {
22: v[j*mat->lda + k] = v[k*mat->lda + j];
23: }
24: }
25: } else {
26: for (k=0;k<n;k++) {
27: for (j=k;j<n;j++) {
28: v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]);
29: }
30: }
31: }
32: MatDenseRestoreArray(A,&v);
33: return 0;
34: }
36: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
37: {
38: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
39: PetscBLASInt info,n;
41: if (!A->rmap->n || !A->cmap->n) return 0;
42: PetscBLASIntCast(A->cmap->n,&n);
43: if (A->factortype == MAT_FACTOR_LU) {
45: if (!mat->fwork) {
46: mat->lfwork = n;
47: PetscMalloc1(mat->lfwork,&mat->fwork);
48: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
49: }
50: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
51: PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
52: PetscFPTrapPop();
53: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
54: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
55: if (A->spd) {
56: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
57: PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
58: PetscFPTrapPop();
59: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
60: #if defined(PETSC_USE_COMPLEX)
61: } else if (A->hermitian) {
64: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
65: PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
66: PetscFPTrapPop();
67: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
68: #endif
69: } else { /* symmetric case */
72: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
73: PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
74: PetscFPTrapPop();
75: MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
76: }
78: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
79: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
81: A->ops->solve = NULL;
82: A->ops->matsolve = NULL;
83: A->ops->solvetranspose = NULL;
84: A->ops->matsolvetranspose = NULL;
85: A->ops->solveadd = NULL;
86: A->ops->solvetransposeadd = NULL;
87: A->factortype = MAT_FACTOR_NONE;
88: PetscFree(A->solvertype);
89: return 0;
90: }
92: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
93: {
94: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
95: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
96: PetscScalar *slot,*bb,*v;
97: const PetscScalar *xx;
99: if (PetscDefined(USE_DEBUG)) {
100: for (i=0; i<N; i++) {
104: }
105: }
106: if (!N) return 0;
108: /* fix right hand side if needed */
109: if (x && b) {
110: Vec xt;
113: VecDuplicate(x,&xt);
114: VecCopy(x,xt);
115: VecScale(xt,-1.0);
116: MatMultAdd(A,xt,b,b);
117: VecDestroy(&xt);
118: VecGetArrayRead(x,&xx);
119: VecGetArray(b,&bb);
120: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
121: VecRestoreArrayRead(x,&xx);
122: VecRestoreArray(b,&bb);
123: }
125: MatDenseGetArray(A,&v);
126: for (i=0; i<N; i++) {
127: slot = v + rows[i]*m;
128: PetscArrayzero(slot,r);
129: }
130: for (i=0; i<N; i++) {
131: slot = v + rows[i];
132: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
133: }
134: if (diag != 0.0) {
136: for (i=0; i<N; i++) {
137: slot = v + (m+1)*rows[i];
138: *slot = diag;
139: }
140: }
141: MatDenseRestoreArray(A,&v);
142: return 0;
143: }
145: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
146: {
147: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
149: if (c->ptapwork) {
150: (*C->ops->matmultnumeric)(A,P,c->ptapwork);
151: (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
152: } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
153: return 0;
154: }
156: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
157: {
158: Mat_SeqDense *c;
159: PetscBool cisdense;
161: MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
162: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
163: if (!cisdense) {
164: PetscBool flg;
166: PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
167: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
168: }
169: MatSetUp(C);
170: c = (Mat_SeqDense*)C->data;
171: MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
172: MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
173: MatSetType(c->ptapwork,((PetscObject)C)->type_name);
174: MatSetUp(c->ptapwork);
175: return 0;
176: }
178: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
179: {
180: Mat B = NULL;
181: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
182: Mat_SeqDense *b;
183: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
184: const MatScalar *av;
185: PetscBool isseqdense;
187: if (reuse == MAT_REUSE_MATRIX) {
188: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
190: }
191: if (reuse != MAT_REUSE_MATRIX) {
192: MatCreate(PetscObjectComm((PetscObject)A),&B);
193: MatSetSizes(B,m,n,m,n);
194: MatSetType(B,MATSEQDENSE);
195: MatSeqDenseSetPreallocation(B,NULL);
196: b = (Mat_SeqDense*)(B->data);
197: } else {
198: b = (Mat_SeqDense*)((*newmat)->data);
199: PetscArrayzero(b->v,m*n);
200: }
201: MatSeqAIJGetArrayRead(A,&av);
202: for (i=0; i<m; i++) {
203: PetscInt j;
204: for (j=0;j<ai[1]-ai[0];j++) {
205: b->v[*aj*m+i] = *av;
206: aj++;
207: av++;
208: }
209: ai++;
210: }
211: MatSeqAIJRestoreArrayRead(A,&av);
213: if (reuse == MAT_INPLACE_MATRIX) {
214: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
216: MatHeaderReplace(A,&B);
217: } else {
218: if (B) *newmat = B;
219: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
221: }
222: return 0;
223: }
225: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
226: {
227: Mat B = NULL;
228: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
229: PetscInt i, j;
230: PetscInt *rows, *nnz;
231: MatScalar *aa = a->v, *vals;
233: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
234: if (reuse != MAT_REUSE_MATRIX) {
235: MatCreate(PetscObjectComm((PetscObject)A),&B);
236: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
237: MatSetType(B,MATSEQAIJ);
238: for (j=0; j<A->cmap->n; j++) {
239: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
240: aa += a->lda;
241: }
242: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
243: } else B = *newmat;
244: aa = a->v;
245: for (j=0; j<A->cmap->n; j++) {
246: PetscInt numRows = 0;
247: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {rows[numRows] = i; vals[numRows++] = aa[i];}
248: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
249: aa += a->lda;
250: }
251: PetscFree3(rows,nnz,vals);
252: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
253: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
255: if (reuse == MAT_INPLACE_MATRIX) {
256: MatHeaderReplace(A,&B);
257: } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
258: return 0;
259: }
261: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
262: {
263: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
264: const PetscScalar *xv;
265: PetscScalar *yv;
266: PetscBLASInt N,m,ldax = 0,lday = 0,one = 1;
268: MatDenseGetArrayRead(X,&xv);
269: MatDenseGetArray(Y,&yv);
270: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
271: PetscBLASIntCast(X->rmap->n,&m);
272: PetscBLASIntCast(x->lda,&ldax);
273: PetscBLASIntCast(y->lda,&lday);
274: if (ldax>m || lday>m) {
275: PetscInt j;
277: for (j=0; j<X->cmap->n; j++) {
278: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
279: }
280: } else {
281: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
282: }
283: MatDenseRestoreArrayRead(X,&xv);
284: MatDenseRestoreArray(Y,&yv);
285: PetscLogFlops(PetscMax(2.0*N-1,0));
286: return 0;
287: }
289: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
290: {
291: PetscLogDouble N = A->rmap->n*A->cmap->n;
293: info->block_size = 1.0;
294: info->nz_allocated = N;
295: info->nz_used = N;
296: info->nz_unneeded = 0;
297: info->assemblies = A->num_ass;
298: info->mallocs = 0;
299: info->memory = ((PetscObject)A)->mem;
300: info->fill_ratio_given = 0;
301: info->fill_ratio_needed = 0;
302: info->factor_mallocs = 0;
303: return 0;
304: }
306: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
307: {
308: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
309: PetscScalar *v;
310: PetscBLASInt one = 1,j,nz,lda = 0;
312: MatDenseGetArray(A,&v);
313: PetscBLASIntCast(a->lda,&lda);
314: if (lda>A->rmap->n) {
315: PetscBLASIntCast(A->rmap->n,&nz);
316: for (j=0; j<A->cmap->n; j++) {
317: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
318: }
319: } else {
320: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
321: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
322: }
323: PetscLogFlops(nz);
324: MatDenseRestoreArray(A,&v);
325: return 0;
326: }
328: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
329: {
330: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
331: PetscInt i,j,m = A->rmap->n,N = a->lda;
332: const PetscScalar *v;
334: *fl = PETSC_FALSE;
335: if (A->rmap->n != A->cmap->n) return 0;
336: MatDenseGetArrayRead(A,&v);
337: for (i=0; i<m; i++) {
338: for (j=i; j<m; j++) {
339: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
340: goto restore;
341: }
342: }
343: }
344: *fl = PETSC_TRUE;
345: restore:
346: MatDenseRestoreArrayRead(A,&v);
347: return 0;
348: }
350: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
351: {
352: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
353: PetscInt i,j,m = A->rmap->n,N = a->lda;
354: const PetscScalar *v;
356: *fl = PETSC_FALSE;
357: if (A->rmap->n != A->cmap->n) return 0;
358: MatDenseGetArrayRead(A,&v);
359: for (i=0; i<m; i++) {
360: for (j=i; j<m; j++) {
361: if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) {
362: goto restore;
363: }
364: }
365: }
366: *fl = PETSC_TRUE;
367: restore:
368: MatDenseRestoreArrayRead(A,&v);
369: return 0;
370: }
372: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
373: {
374: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
375: PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda;
376: PetscBool isdensecpu;
378: PetscLayoutReference(A->rmap,&newi->rmap);
379: PetscLayoutReference(A->cmap,&newi->cmap);
380: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
381: MatDenseSetLDA(newi,lda);
382: }
383: PetscObjectTypeCompare((PetscObject)newi,MATSEQDENSE,&isdensecpu);
384: if (isdensecpu) MatSeqDenseSetPreallocation(newi,NULL);
385: if (cpvalues == MAT_COPY_VALUES) {
386: const PetscScalar *av;
387: PetscScalar *v;
389: MatDenseGetArrayRead(A,&av);
390: MatDenseGetArrayWrite(newi,&v);
391: MatDenseGetLDA(newi,&nlda);
392: m = A->rmap->n;
393: if (lda>m || nlda>m) {
394: for (j=0; j<A->cmap->n; j++) {
395: PetscArraycpy(v+j*nlda,av+j*lda,m);
396: }
397: } else {
398: PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
399: }
400: MatDenseRestoreArrayWrite(newi,&v);
401: MatDenseRestoreArrayRead(A,&av);
402: }
403: return 0;
404: }
406: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
407: {
408: MatCreate(PetscObjectComm((PetscObject)A),newmat);
409: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
410: MatSetType(*newmat,((PetscObject)A)->type_name);
411: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
412: return 0;
413: }
415: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
416: {
417: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
418: PetscBLASInt info;
420: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
421: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(T ? "T" : "N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
422: PetscFPTrapPop();
424: PetscLogFlops(nrhs*(2.0*m*m - m));
425: return 0;
426: }
428: static PetscErrorCode MatConjugate_SeqDense(Mat);
430: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
431: {
432: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
433: PetscBLASInt info;
435: if (A->spd) {
436: if (PetscDefined(USE_COMPLEX) && T) MatConjugate_SeqDense(A);
437: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
438: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
439: PetscFPTrapPop();
441: if (PetscDefined(USE_COMPLEX) && T) MatConjugate_SeqDense(A);
442: #if defined(PETSC_USE_COMPLEX)
443: } else if (A->hermitian) {
444: if (T) MatConjugate_SeqDense(A);
445: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
446: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
447: PetscFPTrapPop();
449: if (T) MatConjugate_SeqDense(A);
450: #endif
451: } else { /* symmetric case */
452: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
453: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
454: PetscFPTrapPop();
456: }
457: PetscLogFlops(nrhs*(2.0*m*m - m));
458: return 0;
459: }
461: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
462: {
463: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
464: PetscBLASInt info;
465: char trans;
467: if (PetscDefined(USE_COMPLEX)) {
468: trans = 'C';
469: } else {
470: trans = 'T';
471: }
472: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
473: PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
474: PetscFPTrapPop();
476: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
477: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
478: PetscFPTrapPop();
480: for (PetscInt j = 0; j < nrhs; j++) {
481: for (PetscInt i = mat->rank; i < k; i++) {
482: x[j*ldx + i] = 0.;
483: }
484: }
485: PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
486: return 0;
487: }
489: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
490: {
491: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
492: PetscBLASInt info;
494: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
495: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
496: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "T", "N", &m,&nrhs,mat->v,&mat->lda,x,&ldx,&info));
497: PetscFPTrapPop();
499: if (PetscDefined(USE_COMPLEX)) MatConjugate_SeqDense(A);
500: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
501: PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info));
502: PetscFPTrapPop();
504: if (PetscDefined(USE_COMPLEX)) MatConjugate_SeqDense(A);
505: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve");
506: PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
507: return 0;
508: }
510: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
511: {
512: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
513: PetscScalar *y;
514: PetscBLASInt m=0, k=0;
516: PetscBLASIntCast(A->rmap->n,&m);
517: PetscBLASIntCast(A->cmap->n,&k);
518: if (k < m) {
519: VecCopy(xx, mat->qrrhs);
520: VecGetArray(mat->qrrhs,&y);
521: } else {
522: VecCopy(xx, yy);
523: VecGetArray(yy,&y);
524: }
525: *_y = y;
526: *_k = k;
527: *_m = m;
528: return 0;
529: }
531: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
532: {
533: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
534: PetscScalar *y = NULL;
535: PetscBLASInt m, k;
537: y = *_y;
538: *_y = NULL;
539: k = *_k;
540: m = *_m;
541: if (k < m) {
542: PetscScalar *yv;
543: VecGetArray(yy,&yv);
544: PetscArraycpy(yv, y, k);
545: VecRestoreArray(yy,&yv);
546: VecRestoreArray(mat->qrrhs, &y);
547: } else {
548: VecRestoreArray(yy,&y);
549: }
550: return 0;
551: }
553: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
554: {
555: PetscScalar *y = NULL;
556: PetscBLASInt m = 0, k = 0;
558: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
559: MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE);
560: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
561: return 0;
562: }
564: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
565: {
566: PetscScalar *y = NULL;
567: PetscBLASInt m = 0, k = 0;
569: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
570: MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE);
571: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
572: return 0;
573: }
575: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
576: {
577: PetscScalar *y = NULL;
578: PetscBLASInt m = 0, k = 0;
580: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
581: MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE);
582: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
583: return 0;
584: }
586: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
587: {
588: PetscScalar *y = NULL;
589: PetscBLASInt m = 0, k = 0;
591: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
592: MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE);
593: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
594: return 0;
595: }
597: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
598: {
599: PetscScalar *y = NULL;
600: PetscBLASInt m = 0, k = 0;
602: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
603: MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
604: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
605: return 0;
606: }
608: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
609: {
610: PetscScalar *y = NULL;
611: PetscBLASInt m = 0, k = 0;
613: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
614: MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
615: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
616: return 0;
617: }
619: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
620: {
621: const PetscScalar *b;
622: PetscScalar *y;
623: PetscInt n, _ldb, _ldx;
624: PetscBLASInt nrhs=0,m=0,k=0,ldb=0,ldx=0,ldy=0;
626: *_ldy=0; *_m=0; *_nrhs=0; *_k=0; *_y = NULL;
627: PetscBLASIntCast(A->rmap->n,&m);
628: PetscBLASIntCast(A->cmap->n,&k);
629: MatGetSize(B,NULL,&n);
630: PetscBLASIntCast(n,&nrhs);
631: MatDenseGetLDA(B,&_ldb);
632: PetscBLASIntCast(_ldb, &ldb);
633: MatDenseGetLDA(X,&_ldx);
634: PetscBLASIntCast(_ldx, &ldx);
635: if (ldx < m) {
636: MatDenseGetArrayRead(B,&b);
637: PetscMalloc1(nrhs * m, &y);
638: if (ldb == m) {
639: PetscArraycpy(y,b,ldb*nrhs);
640: } else {
641: for (PetscInt j = 0; j < nrhs; j++) {
642: PetscArraycpy(&y[j*m],&b[j*ldb],m);
643: }
644: }
645: ldy = m;
646: MatDenseRestoreArrayRead(B,&b);
647: } else {
648: if (ldb == ldx) {
649: MatCopy(B, X, SAME_NONZERO_PATTERN);
650: MatDenseGetArray(X,&y);
651: } else {
652: MatDenseGetArray(X,&y);
653: MatDenseGetArrayRead(B,&b);
654: for (PetscInt j = 0; j < nrhs; j++) {
655: PetscArraycpy(&y[j*ldx],&b[j*ldb],m);
656: }
657: MatDenseRestoreArrayRead(B,&b);
658: }
659: ldy = ldx;
660: }
661: *_y = y;
662: *_ldy = ldy;
663: *_k = k;
664: *_m = m;
665: *_nrhs = nrhs;
666: return 0;
667: }
669: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
670: {
671: PetscScalar *y;
672: PetscInt _ldx;
673: PetscBLASInt k,ldy,nrhs,ldx=0;
675: y = *_y;
676: *_y = NULL;
677: k = *_k;
678: ldy = *_ldy;
679: nrhs = *_nrhs;
680: MatDenseGetLDA(X,&_ldx);
681: PetscBLASIntCast(_ldx, &ldx);
682: if (ldx != ldy) {
683: PetscScalar *xv;
684: MatDenseGetArray(X,&xv);
685: for (PetscInt j = 0; j < nrhs; j++) {
686: PetscArraycpy(&xv[j*ldx],&y[j*ldy],k);
687: }
688: MatDenseRestoreArray(X,&xv);
689: PetscFree(y);
690: } else {
691: MatDenseRestoreArray(X,&y);
692: }
693: return 0;
694: }
696: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
697: {
698: PetscScalar *y;
699: PetscBLASInt m, k, ldy, nrhs;
701: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
702: MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE);
703: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
704: return 0;
705: }
707: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
708: {
709: PetscScalar *y;
710: PetscBLASInt m, k, ldy, nrhs;
712: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
713: MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE);
714: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
715: return 0;
716: }
718: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
719: {
720: PetscScalar *y;
721: PetscBLASInt m, k, ldy, nrhs;
723: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
724: MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE);
725: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
726: return 0;
727: }
729: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
730: {
731: PetscScalar *y;
732: PetscBLASInt m, k, ldy, nrhs;
734: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
735: MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE);
736: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
737: return 0;
738: }
740: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
741: {
742: PetscScalar *y;
743: PetscBLASInt m, k, ldy, nrhs;
745: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
746: MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
747: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
748: return 0;
749: }
751: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
752: {
753: PetscScalar *y;
754: PetscBLASInt m, k, ldy, nrhs;
756: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
757: MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
758: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k);
759: return 0;
760: }
762: static PetscErrorCode MatConjugate_SeqDense(Mat);
764: /* ---------------------------------------------------------------*/
765: /* COMMENT: I have chosen to hide row permutation in the pivots,
766: rather than put it in the Mat->row slot.*/
767: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
768: {
769: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
770: PetscBLASInt n,m,info;
772: PetscBLASIntCast(A->cmap->n,&n);
773: PetscBLASIntCast(A->rmap->n,&m);
774: if (!mat->pivots) {
775: PetscMalloc1(A->rmap->n,&mat->pivots);
776: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
777: }
778: if (!A->rmap->n || !A->cmap->n) return 0;
779: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
780: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
781: PetscFPTrapPop();
786: A->ops->solve = MatSolve_SeqDense_LU;
787: A->ops->matsolve = MatMatSolve_SeqDense_LU;
788: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
789: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
790: A->factortype = MAT_FACTOR_LU;
792: PetscFree(A->solvertype);
793: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
795: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
796: return 0;
797: }
799: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
800: {
801: MatFactorInfo info;
803: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
804: (*fact->ops->lufactor)(fact,NULL,NULL,&info);
805: return 0;
806: }
808: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
809: {
810: fact->preallocated = PETSC_TRUE;
811: fact->assembled = PETSC_TRUE;
812: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
813: return 0;
814: }
816: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
817: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
818: {
819: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
820: PetscBLASInt info,n;
822: PetscBLASIntCast(A->cmap->n,&n);
823: if (!A->rmap->n || !A->cmap->n) return 0;
824: if (A->spd) {
825: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
826: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
827: PetscFPTrapPop();
828: #if defined(PETSC_USE_COMPLEX)
829: } else if (A->hermitian) {
830: if (!mat->pivots) {
831: PetscMalloc1(A->rmap->n,&mat->pivots);
832: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
833: }
834: if (!mat->fwork) {
835: PetscScalar dummy;
837: mat->lfwork = -1;
838: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
839: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
840: PetscFPTrapPop();
841: mat->lfwork = (PetscInt)PetscRealPart(dummy);
842: PetscMalloc1(mat->lfwork,&mat->fwork);
843: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
844: }
845: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
846: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
847: PetscFPTrapPop();
848: #endif
849: } else { /* symmetric case */
850: if (!mat->pivots) {
851: PetscMalloc1(A->rmap->n,&mat->pivots);
852: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
853: }
854: if (!mat->fwork) {
855: PetscScalar dummy;
857: mat->lfwork = -1;
858: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
859: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
860: PetscFPTrapPop();
861: mat->lfwork = (PetscInt)PetscRealPart(dummy);
862: PetscMalloc1(mat->lfwork,&mat->fwork);
863: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
864: }
865: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
866: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
867: PetscFPTrapPop();
868: }
871: A->ops->solve = MatSolve_SeqDense_Cholesky;
872: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
873: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
874: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
875: A->factortype = MAT_FACTOR_CHOLESKY;
877: PetscFree(A->solvertype);
878: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
880: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
881: return 0;
882: }
884: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
885: {
886: MatFactorInfo info;
888: info.fill = 1.0;
890: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
891: (*fact->ops->choleskyfactor)(fact,NULL,&info);
892: return 0;
893: }
895: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
896: {
897: fact->assembled = PETSC_TRUE;
898: fact->preallocated = PETSC_TRUE;
899: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
900: return 0;
901: }
903: PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
904: {
905: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
906: PetscBLASInt n,m,info, min, max;
908: PetscBLASIntCast(A->cmap->n,&n);
909: PetscBLASIntCast(A->rmap->n,&m);
910: max = PetscMax(m, n);
911: min = PetscMin(m, n);
912: if (!mat->tau) {
913: PetscMalloc1(min,&mat->tau);
914: PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
915: }
916: if (!mat->pivots) {
917: PetscMalloc1(n,&mat->pivots);
918: PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar));
919: }
920: if (!mat->qrrhs) {
921: MatCreateVecs(A, NULL, &(mat->qrrhs));
922: }
923: if (!A->rmap->n || !A->cmap->n) return 0;
924: if (!mat->fwork) {
925: PetscScalar dummy;
927: mat->lfwork = -1;
928: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
929: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
930: PetscFPTrapPop();
931: mat->lfwork = (PetscInt)PetscRealPart(dummy);
932: PetscMalloc1(mat->lfwork,&mat->fwork);
933: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
934: }
935: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
936: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
937: PetscFPTrapPop();
939: // TODO: try to estimate rank or test for and use geqp3 for rank revealing QR. For now just say rank is min of m and n
940: mat->rank = min;
942: A->ops->solve = MatSolve_SeqDense_QR;
943: A->ops->matsolve = MatMatSolve_SeqDense_QR;
944: A->factortype = MAT_FACTOR_QR;
945: if (m == n) {
946: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
947: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
948: }
950: PetscFree(A->solvertype);
951: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
953: PetscLogFlops(2.0*min*min*(max-min/3.0));
954: return 0;
955: }
957: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
958: {
959: MatFactorInfo info;
961: info.fill = 1.0;
963: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
964: PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info));
965: return 0;
966: }
968: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
969: {
970: fact->assembled = PETSC_TRUE;
971: fact->preallocated = PETSC_TRUE;
972: PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
973: return 0;
974: }
976: /* uses LAPACK */
977: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
978: {
979: MatCreate(PetscObjectComm((PetscObject)A),fact);
980: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
981: MatSetType(*fact,MATDENSE);
982: (*fact)->trivialsymbolic = PETSC_TRUE;
983: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
984: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
985: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
986: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
987: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
988: } else if (ftype == MAT_FACTOR_QR) {
989: PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
990: }
991: (*fact)->factortype = ftype;
993: PetscFree((*fact)->solvertype);
994: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
995: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU]);
996: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU]);
997: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]);
998: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC]);
999: return 0;
1000: }
1002: /* ------------------------------------------------------------------*/
1003: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1004: {
1005: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1006: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
1007: const PetscScalar *b;
1008: PetscInt m = A->rmap->n,i;
1009: PetscBLASInt o = 1,bm = 0;
1011: #if defined(PETSC_HAVE_CUDA)
1013: #endif
1014: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1015: PetscBLASIntCast(m,&bm);
1016: if (flag & SOR_ZERO_INITIAL_GUESS) {
1017: /* this is a hack fix, should have another version without the second BLASdotu */
1018: VecSet(xx,zero);
1019: }
1020: VecGetArray(xx,&x);
1021: VecGetArrayRead(bb,&b);
1022: its = its*lits;
1024: while (its--) {
1025: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1026: for (i=0; i<m; i++) {
1027: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1028: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1029: }
1030: }
1031: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1032: for (i=m-1; i>=0; i--) {
1033: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1034: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1035: }
1036: }
1037: }
1038: VecRestoreArrayRead(bb,&b);
1039: VecRestoreArray(xx,&x);
1040: return 0;
1041: }
1043: /* -----------------------------------------------------------------*/
1044: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1045: {
1046: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1047: const PetscScalar *v = mat->v,*x;
1048: PetscScalar *y;
1049: PetscBLASInt m, n,_One=1;
1050: PetscScalar _DOne=1.0,_DZero=0.0;
1052: PetscBLASIntCast(A->rmap->n,&m);
1053: PetscBLASIntCast(A->cmap->n,&n);
1054: VecGetArrayRead(xx,&x);
1055: VecGetArrayWrite(yy,&y);
1056: if (!A->rmap->n || !A->cmap->n) {
1057: PetscBLASInt i;
1058: for (i=0; i<n; i++) y[i] = 0.0;
1059: } else {
1060: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1061: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1062: }
1063: VecRestoreArrayRead(xx,&x);
1064: VecRestoreArrayWrite(yy,&y);
1065: return 0;
1066: }
1068: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1069: {
1070: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1071: PetscScalar *y,_DOne=1.0,_DZero=0.0;
1072: PetscBLASInt m, n, _One=1;
1073: const PetscScalar *v = mat->v,*x;
1075: PetscBLASIntCast(A->rmap->n,&m);
1076: PetscBLASIntCast(A->cmap->n,&n);
1077: VecGetArrayRead(xx,&x);
1078: VecGetArrayWrite(yy,&y);
1079: if (!A->rmap->n || !A->cmap->n) {
1080: PetscBLASInt i;
1081: for (i=0; i<m; i++) y[i] = 0.0;
1082: } else {
1083: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1084: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1085: }
1086: VecRestoreArrayRead(xx,&x);
1087: VecRestoreArrayWrite(yy,&y);
1088: return 0;
1089: }
1091: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1092: {
1093: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1094: const PetscScalar *v = mat->v,*x;
1095: PetscScalar *y,_DOne=1.0;
1096: PetscBLASInt m, n, _One=1;
1098: PetscBLASIntCast(A->rmap->n,&m);
1099: PetscBLASIntCast(A->cmap->n,&n);
1100: VecCopy(zz,yy);
1101: if (!A->rmap->n || !A->cmap->n) return 0;
1102: VecGetArrayRead(xx,&x);
1103: VecGetArray(yy,&y);
1104: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1105: VecRestoreArrayRead(xx,&x);
1106: VecRestoreArray(yy,&y);
1107: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1108: return 0;
1109: }
1111: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1112: {
1113: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1114: const PetscScalar *v = mat->v,*x;
1115: PetscScalar *y;
1116: PetscBLASInt m, n, _One=1;
1117: PetscScalar _DOne=1.0;
1119: PetscBLASIntCast(A->rmap->n,&m);
1120: PetscBLASIntCast(A->cmap->n,&n);
1121: VecCopy(zz,yy);
1122: if (!A->rmap->n || !A->cmap->n) return 0;
1123: VecGetArrayRead(xx,&x);
1124: VecGetArray(yy,&y);
1125: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1126: VecRestoreArrayRead(xx,&x);
1127: VecRestoreArray(yy,&y);
1128: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1129: return 0;
1130: }
1132: /* -----------------------------------------------------------------*/
1133: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1134: {
1135: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1136: PetscInt i;
1138: *ncols = A->cmap->n;
1139: if (cols) {
1140: PetscMalloc1(A->cmap->n,cols);
1141: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1142: }
1143: if (vals) {
1144: const PetscScalar *v;
1146: MatDenseGetArrayRead(A,&v);
1147: PetscMalloc1(A->cmap->n,vals);
1148: v += row;
1149: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1150: MatDenseRestoreArrayRead(A,&v);
1151: }
1152: return 0;
1153: }
1155: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1156: {
1157: if (ncols) *ncols = 0;
1158: if (cols) PetscFree(*cols);
1159: if (vals) PetscFree(*vals);
1160: return 0;
1161: }
1162: /* ----------------------------------------------------------------*/
1163: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1164: {
1165: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1166: PetscScalar *av;
1167: PetscInt i,j,idx=0;
1168: #if defined(PETSC_HAVE_CUDA)
1169: PetscOffloadMask oldf;
1170: #endif
1172: MatDenseGetArray(A,&av);
1173: if (!mat->roworiented) {
1174: if (addv == INSERT_VALUES) {
1175: for (j=0; j<n; j++) {
1176: if (indexn[j] < 0) {idx += m; continue;}
1178: for (i=0; i<m; i++) {
1179: if (indexm[i] < 0) {idx++; continue;}
1181: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1182: }
1183: }
1184: } else {
1185: for (j=0; j<n; j++) {
1186: if (indexn[j] < 0) {idx += m; continue;}
1188: for (i=0; i<m; i++) {
1189: if (indexm[i] < 0) {idx++; continue;}
1191: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1192: }
1193: }
1194: }
1195: } else {
1196: if (addv == INSERT_VALUES) {
1197: for (i=0; i<m; i++) {
1198: if (indexm[i] < 0) { idx += n; continue;}
1200: for (j=0; j<n; j++) {
1201: if (indexn[j] < 0) { idx++; continue;}
1203: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1204: }
1205: }
1206: } else {
1207: for (i=0; i<m; i++) {
1208: if (indexm[i] < 0) { idx += n; continue;}
1210: for (j=0; j<n; j++) {
1211: if (indexn[j] < 0) { idx++; continue;}
1213: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1214: }
1215: }
1216: }
1217: }
1218: /* hack to prevent unneeded copy to the GPU while returning the array */
1219: #if defined(PETSC_HAVE_CUDA)
1220: oldf = A->offloadmask;
1221: A->offloadmask = PETSC_OFFLOAD_GPU;
1222: #endif
1223: MatDenseRestoreArray(A,&av);
1224: #if defined(PETSC_HAVE_CUDA)
1225: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1226: #endif
1227: return 0;
1228: }
1230: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1231: {
1232: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1233: const PetscScalar *vv;
1234: PetscInt i,j;
1236: MatDenseGetArrayRead(A,&vv);
1237: /* row-oriented output */
1238: for (i=0; i<m; i++) {
1239: if (indexm[i] < 0) {v += n;continue;}
1241: for (j=0; j<n; j++) {
1242: if (indexn[j] < 0) {v++; continue;}
1244: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1245: }
1246: }
1247: MatDenseRestoreArrayRead(A,&vv);
1248: return 0;
1249: }
1251: /* -----------------------------------------------------------------*/
1253: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1254: {
1255: PetscBool skipHeader;
1256: PetscViewerFormat format;
1257: PetscInt header[4],M,N,m,lda,i,j,k;
1258: const PetscScalar *v;
1259: PetscScalar *vwork;
1261: PetscViewerSetUp(viewer);
1262: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1263: PetscViewerGetFormat(viewer,&format);
1264: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1266: MatGetSize(mat,&M,&N);
1268: /* write matrix header */
1269: header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1270: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1271: if (!skipHeader) PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);
1273: MatGetLocalSize(mat,&m,NULL);
1274: if (format != PETSC_VIEWER_NATIVE) {
1275: PetscInt nnz = m*N, *iwork;
1276: /* store row lengths for each row */
1277: PetscMalloc1(nnz,&iwork);
1278: for (i=0; i<m; i++) iwork[i] = N;
1279: PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1280: /* store column indices (zero start index) */
1281: for (k=0, i=0; i<m; i++)
1282: for (j=0; j<N; j++, k++)
1283: iwork[k] = j;
1284: PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1285: PetscFree(iwork);
1286: }
1287: /* store matrix values as a dense matrix in row major order */
1288: PetscMalloc1(m*N,&vwork);
1289: MatDenseGetArrayRead(mat,&v);
1290: MatDenseGetLDA(mat,&lda);
1291: for (k=0, i=0; i<m; i++)
1292: for (j=0; j<N; j++, k++)
1293: vwork[k] = v[i+lda*j];
1294: MatDenseRestoreArrayRead(mat,&v);
1295: PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1296: PetscFree(vwork);
1297: return 0;
1298: }
1300: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1301: {
1302: PetscBool skipHeader;
1303: PetscInt header[4],M,N,m,nz,lda,i,j,k;
1304: PetscInt rows,cols;
1305: PetscScalar *v,*vwork;
1307: PetscViewerSetUp(viewer);
1308: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1310: if (!skipHeader) {
1311: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1313: M = header[1]; N = header[2];
1316: nz = header[3];
1318: } else {
1319: MatGetSize(mat,&M,&N);
1321: nz = MATRIX_BINARY_FORMAT_DENSE;
1322: }
1324: /* setup global sizes if not set */
1325: if (mat->rmap->N < 0) mat->rmap->N = M;
1326: if (mat->cmap->N < 0) mat->cmap->N = N;
1327: MatSetUp(mat);
1328: /* check if global sizes are correct */
1329: MatGetSize(mat,&rows,&cols);
1332: MatGetSize(mat,NULL,&N);
1333: MatGetLocalSize(mat,&m,NULL);
1334: MatDenseGetArray(mat,&v);
1335: MatDenseGetLDA(mat,&lda);
1336: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1337: PetscInt nnz = m*N;
1338: /* read in matrix values */
1339: PetscMalloc1(nnz,&vwork);
1340: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1341: /* store values in column major order */
1342: for (j=0; j<N; j++)
1343: for (i=0; i<m; i++)
1344: v[i+lda*j] = vwork[i*N+j];
1345: PetscFree(vwork);
1346: } else { /* matrix in file is sparse format */
1347: PetscInt nnz = 0, *rlens, *icols;
1348: /* read in row lengths */
1349: PetscMalloc1(m,&rlens);
1350: PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1351: for (i=0; i<m; i++) nnz += rlens[i];
1352: /* read in column indices and values */
1353: PetscMalloc2(nnz,&icols,nnz,&vwork);
1354: PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1355: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1356: /* store values in column major order */
1357: for (k=0, i=0; i<m; i++)
1358: for (j=0; j<rlens[i]; j++, k++)
1359: v[i+lda*icols[k]] = vwork[k];
1360: PetscFree(rlens);
1361: PetscFree2(icols,vwork);
1362: }
1363: MatDenseRestoreArray(mat,&v);
1364: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1365: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1366: return 0;
1367: }
1369: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1370: {
1371: PetscBool isbinary, ishdf5;
1375: /* force binary viewer to load .info file if it has not yet done so */
1376: PetscViewerSetUp(viewer);
1377: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1378: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1379: if (isbinary) {
1380: MatLoad_Dense_Binary(newMat,viewer);
1381: } else if (ishdf5) {
1382: #if defined(PETSC_HAVE_HDF5)
1383: MatLoad_Dense_HDF5(newMat,viewer);
1384: #else
1385: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1386: #endif
1387: } else {
1388: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1389: }
1390: return 0;
1391: }
1393: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1394: {
1395: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1396: PetscInt i,j;
1397: const char *name;
1398: PetscScalar *v,*av;
1399: PetscViewerFormat format;
1400: #if defined(PETSC_USE_COMPLEX)
1401: PetscBool allreal = PETSC_TRUE;
1402: #endif
1404: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1405: PetscViewerGetFormat(viewer,&format);
1406: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1407: return 0; /* do nothing for now */
1408: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1409: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1410: for (i=0; i<A->rmap->n; i++) {
1411: v = av + i;
1412: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
1413: for (j=0; j<A->cmap->n; j++) {
1414: #if defined(PETSC_USE_COMPLEX)
1415: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1416: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1417: } else if (PetscRealPart(*v)) {
1418: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)PetscRealPart(*v));
1419: }
1420: #else
1421: if (*v) {
1422: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)*v);
1423: }
1424: #endif
1425: v += a->lda;
1426: }
1427: PetscViewerASCIIPrintf(viewer,"\n");
1428: }
1429: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1430: } else {
1431: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1432: #if defined(PETSC_USE_COMPLEX)
1433: /* determine if matrix has all real values */
1434: v = av;
1435: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1436: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1437: }
1438: #endif
1439: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1440: PetscObjectGetName((PetscObject)A,&name);
1441: PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",A->rmap->n,A->cmap->n);
1442: PetscViewerASCIIPrintf(viewer,"%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n",name,A->rmap->n,A->cmap->n);
1443: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1444: }
1446: for (i=0; i<A->rmap->n; i++) {
1447: v = av + i;
1448: for (j=0; j<A->cmap->n; j++) {
1449: #if defined(PETSC_USE_COMPLEX)
1450: if (allreal) {
1451: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1452: } else {
1453: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1454: }
1455: #else
1456: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1457: #endif
1458: v += a->lda;
1459: }
1460: PetscViewerASCIIPrintf(viewer,"\n");
1461: }
1462: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1463: PetscViewerASCIIPrintf(viewer,"];\n");
1464: }
1465: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1466: }
1467: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1468: PetscViewerFlush(viewer);
1469: return 0;
1470: }
1472: #include <petscdraw.h>
1473: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1474: {
1475: Mat A = (Mat) Aa;
1476: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1477: int color = PETSC_DRAW_WHITE;
1478: const PetscScalar *v;
1479: PetscViewer viewer;
1480: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1481: PetscViewerFormat format;
1482: PetscErrorCode ierr;
1484: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1485: PetscViewerGetFormat(viewer,&format);
1486: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1488: /* Loop over matrix elements drawing boxes */
1489: MatDenseGetArrayRead(A,&v);
1490: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1491: PetscDrawCollectiveBegin(draw);
1492: /* Blue for negative and Red for positive */
1493: for (j = 0; j < n; j++) {
1494: x_l = j; x_r = x_l + 1.0;
1495: for (i = 0; i < m; i++) {
1496: y_l = m - i - 1.0;
1497: y_r = y_l + 1.0;
1498: if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED;
1499: else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE;
1500: else continue;
1501: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1502: }
1503: }
1504: PetscDrawCollectiveEnd(draw);
1505: } else {
1506: /* use contour shading to indicate magnitude of values */
1507: /* first determine max of all nonzero values */
1508: PetscReal minv = 0.0, maxv = 0.0;
1509: PetscDraw popup;
1511: for (i=0; i < m*n; i++) {
1512: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1513: }
1514: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1515: PetscDrawGetPopup(draw,&popup);
1516: PetscDrawScalePopup(popup,minv,maxv);
1518: PetscDrawCollectiveBegin(draw);
1519: for (j=0; j<n; j++) {
1520: x_l = j;
1521: x_r = x_l + 1.0;
1522: for (i=0; i<m; i++) {
1523: y_l = m - i - 1.0;
1524: y_r = y_l + 1.0;
1525: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1526: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1527: }
1528: }
1529: PetscDrawCollectiveEnd(draw);
1530: }
1531: MatDenseRestoreArrayRead(A,&v);
1532: return 0;
1533: }
1535: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1536: {
1537: PetscDraw draw;
1538: PetscBool isnull;
1539: PetscReal xr,yr,xl,yl,h,w;
1541: PetscViewerDrawGetDraw(viewer,0,&draw);
1542: PetscDrawIsNull(draw,&isnull);
1543: if (isnull) return 0;
1545: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1546: xr += w; yr += h; xl = -w; yl = -h;
1547: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1548: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1549: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1550: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1551: PetscDrawSave(draw);
1552: return 0;
1553: }
1555: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1556: {
1557: PetscBool iascii,isbinary,isdraw;
1559: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1560: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1561: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1562: if (iascii) {
1563: MatView_SeqDense_ASCII(A,viewer);
1564: } else if (isbinary) {
1565: MatView_Dense_Binary(A,viewer);
1566: } else if (isdraw) {
1567: MatView_SeqDense_Draw(A,viewer);
1568: }
1569: return 0;
1570: }
1572: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1573: {
1574: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1579: a->unplacedarray = a->v;
1580: a->unplaced_user_alloc = a->user_alloc;
1581: a->v = (PetscScalar*) array;
1582: a->user_alloc = PETSC_TRUE;
1583: #if defined(PETSC_HAVE_CUDA)
1584: A->offloadmask = PETSC_OFFLOAD_CPU;
1585: #endif
1586: return 0;
1587: }
1589: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1590: {
1591: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1595: a->v = a->unplacedarray;
1596: a->user_alloc = a->unplaced_user_alloc;
1597: a->unplacedarray = NULL;
1598: #if defined(PETSC_HAVE_CUDA)
1599: A->offloadmask = PETSC_OFFLOAD_CPU;
1600: #endif
1601: return 0;
1602: }
1604: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1605: {
1606: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1610: if (!a->user_alloc) PetscFree(a->v);
1611: a->v = (PetscScalar*) array;
1612: a->user_alloc = PETSC_FALSE;
1613: #if defined(PETSC_HAVE_CUDA)
1614: A->offloadmask = PETSC_OFFLOAD_CPU;
1615: #endif
1616: return 0;
1617: }
1619: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1620: {
1621: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1623: #if defined(PETSC_USE_LOG)
1624: PetscLogObjectState((PetscObject)mat,"Rows %" PetscInt_FMT " Cols %" PetscInt_FMT,mat->rmap->n,mat->cmap->n);
1625: #endif
1626: VecDestroy(&(l->qrrhs));
1627: PetscFree(l->tau);
1628: PetscFree(l->pivots);
1629: PetscFree(l->fwork);
1630: MatDestroy(&l->ptapwork);
1631: if (!l->user_alloc) PetscFree(l->v);
1632: if (!l->unplaced_user_alloc) PetscFree(l->unplacedarray);
1635: VecDestroy(&l->cvec);
1636: MatDestroy(&l->cmat);
1637: PetscFree(mat->data);
1639: PetscObjectChangeTypeName((PetscObject)mat,NULL);
1640: PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);
1641: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1642: PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1643: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1644: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1645: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1646: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1647: PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1648: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1649: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1650: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1651: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1652: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1653: #if defined(PETSC_HAVE_ELEMENTAL)
1654: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1655: #endif
1656: #if defined(PETSC_HAVE_SCALAPACK)
1657: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1658: #endif
1659: #if defined(PETSC_HAVE_CUDA)
1660: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1661: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1662: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1663: #endif
1664: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1665: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1666: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1667: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1668: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);
1670: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1671: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1672: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1673: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1674: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1675: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1676: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1677: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1678: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1679: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1680: return 0;
1681: }
1683: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1684: {
1685: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1686: PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1687: PetscScalar *v,tmp;
1689: if (reuse == MAT_INPLACE_MATRIX) {
1690: if (m == n) { /* in place transpose */
1691: MatDenseGetArray(A,&v);
1692: for (j=0; j<m; j++) {
1693: for (k=0; k<j; k++) {
1694: tmp = v[j + k*M];
1695: v[j + k*M] = v[k + j*M];
1696: v[k + j*M] = tmp;
1697: }
1698: }
1699: MatDenseRestoreArray(A,&v);
1700: } else { /* reuse memory, temporary allocates new memory */
1701: PetscScalar *v2;
1702: PetscLayout tmplayout;
1704: PetscMalloc1((size_t)m*n,&v2);
1705: MatDenseGetArray(A,&v);
1706: for (j=0; j<n; j++) {
1707: for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1708: }
1709: PetscArraycpy(v,v2,(size_t)m*n);
1710: PetscFree(v2);
1711: MatDenseRestoreArray(A,&v);
1712: /* cleanup size dependent quantities */
1713: VecDestroy(&mat->cvec);
1714: MatDestroy(&mat->cmat);
1715: PetscFree(mat->pivots);
1716: PetscFree(mat->fwork);
1717: MatDestroy(&mat->ptapwork);
1718: /* swap row/col layouts */
1719: mat->lda = n;
1720: tmplayout = A->rmap;
1721: A->rmap = A->cmap;
1722: A->cmap = tmplayout;
1723: }
1724: } else { /* out-of-place transpose */
1725: Mat tmat;
1726: Mat_SeqDense *tmatd;
1727: PetscScalar *v2;
1728: PetscInt M2;
1730: if (reuse == MAT_INITIAL_MATRIX) {
1731: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1732: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1733: MatSetType(tmat,((PetscObject)A)->type_name);
1734: MatSeqDenseSetPreallocation(tmat,NULL);
1735: } else tmat = *matout;
1737: MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1738: MatDenseGetArray(tmat,&v2);
1739: tmatd = (Mat_SeqDense*)tmat->data;
1740: M2 = tmatd->lda;
1741: for (j=0; j<n; j++) {
1742: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1743: }
1744: MatDenseRestoreArray(tmat,&v2);
1745: MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1746: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1747: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1748: *matout = tmat;
1749: }
1750: return 0;
1751: }
1753: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1754: {
1755: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1756: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1757: PetscInt i;
1758: const PetscScalar *v1,*v2;
1760: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return 0;}
1761: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return 0;}
1762: MatDenseGetArrayRead(A1,&v1);
1763: MatDenseGetArrayRead(A2,&v2);
1764: for (i=0; i<A1->cmap->n; i++) {
1765: PetscArraycmp(v1,v2,A1->rmap->n,flg);
1766: if (*flg == PETSC_FALSE) return 0;
1767: v1 += mat1->lda;
1768: v2 += mat2->lda;
1769: }
1770: MatDenseRestoreArrayRead(A1,&v1);
1771: MatDenseRestoreArrayRead(A2,&v2);
1772: *flg = PETSC_TRUE;
1773: return 0;
1774: }
1776: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1777: {
1778: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1779: PetscInt i,n,len;
1780: PetscScalar *x;
1781: const PetscScalar *vv;
1783: VecGetSize(v,&n);
1784: VecGetArray(v,&x);
1785: len = PetscMin(A->rmap->n,A->cmap->n);
1786: MatDenseGetArrayRead(A,&vv);
1788: for (i=0; i<len; i++) {
1789: x[i] = vv[i*mat->lda + i];
1790: }
1791: MatDenseRestoreArrayRead(A,&vv);
1792: VecRestoreArray(v,&x);
1793: return 0;
1794: }
1796: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1797: {
1798: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1799: const PetscScalar *l,*r;
1800: PetscScalar x,*v,*vv;
1801: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1803: MatDenseGetArray(A,&vv);
1804: if (ll) {
1805: VecGetSize(ll,&m);
1806: VecGetArrayRead(ll,&l);
1808: for (i=0; i<m; i++) {
1809: x = l[i];
1810: v = vv + i;
1811: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1812: }
1813: VecRestoreArrayRead(ll,&l);
1814: PetscLogFlops(1.0*n*m);
1815: }
1816: if (rr) {
1817: VecGetSize(rr,&n);
1818: VecGetArrayRead(rr,&r);
1820: for (i=0; i<n; i++) {
1821: x = r[i];
1822: v = vv + i*mat->lda;
1823: for (j=0; j<m; j++) (*v++) *= x;
1824: }
1825: VecRestoreArrayRead(rr,&r);
1826: PetscLogFlops(1.0*n*m);
1827: }
1828: MatDenseRestoreArray(A,&vv);
1829: return 0;
1830: }
1832: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1833: {
1834: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1835: PetscScalar *v,*vv;
1836: PetscReal sum = 0.0;
1837: PetscInt lda, m=A->rmap->n,i,j;
1839: MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1840: MatDenseGetLDA(A,&lda);
1841: v = vv;
1842: if (type == NORM_FROBENIUS) {
1843: if (lda>m) {
1844: for (j=0; j<A->cmap->n; j++) {
1845: v = vv+j*lda;
1846: for (i=0; i<m; i++) {
1847: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1848: }
1849: }
1850: } else {
1851: #if defined(PETSC_USE_REAL___FP16)
1852: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1853: PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one));
1854: }
1855: #else
1856: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1857: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1858: }
1859: }
1860: *nrm = PetscSqrtReal(sum);
1861: #endif
1862: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1863: } else if (type == NORM_1) {
1864: *nrm = 0.0;
1865: for (j=0; j<A->cmap->n; j++) {
1866: v = vv + j*mat->lda;
1867: sum = 0.0;
1868: for (i=0; i<A->rmap->n; i++) {
1869: sum += PetscAbsScalar(*v); v++;
1870: }
1871: if (sum > *nrm) *nrm = sum;
1872: }
1873: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1874: } else if (type == NORM_INFINITY) {
1875: *nrm = 0.0;
1876: for (j=0; j<A->rmap->n; j++) {
1877: v = vv + j;
1878: sum = 0.0;
1879: for (i=0; i<A->cmap->n; i++) {
1880: sum += PetscAbsScalar(*v); v += mat->lda;
1881: }
1882: if (sum > *nrm) *nrm = sum;
1883: }
1884: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1885: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1886: MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1887: return 0;
1888: }
1890: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1891: {
1892: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1894: switch (op) {
1895: case MAT_ROW_ORIENTED:
1896: aij->roworiented = flg;
1897: break;
1898: case MAT_NEW_NONZERO_LOCATIONS:
1899: case MAT_NEW_NONZERO_LOCATION_ERR:
1900: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1901: case MAT_FORCE_DIAGONAL_ENTRIES:
1902: case MAT_KEEP_NONZERO_PATTERN:
1903: case MAT_IGNORE_OFF_PROC_ENTRIES:
1904: case MAT_USE_HASH_TABLE:
1905: case MAT_IGNORE_ZERO_ENTRIES:
1906: case MAT_IGNORE_LOWER_TRIANGULAR:
1907: case MAT_SORTED_FULL:
1908: PetscInfo(A,"Option %s ignored\n",MatOptions[op]);
1909: break;
1910: case MAT_SPD:
1911: case MAT_SYMMETRIC:
1912: case MAT_STRUCTURALLY_SYMMETRIC:
1913: case MAT_HERMITIAN:
1914: case MAT_SYMMETRY_ETERNAL:
1915: /* These options are handled directly by MatSetOption() */
1916: break;
1917: default:
1918: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1919: }
1920: return 0;
1921: }
1923: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1924: {
1925: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1926: PetscInt lda=l->lda,m=A->rmap->n,n=A->cmap->n,j;
1927: PetscScalar *v;
1929: MatDenseGetArrayWrite(A,&v);
1930: if (lda>m) {
1931: for (j=0; j<n; j++) {
1932: PetscArrayzero(v+j*lda,m);
1933: }
1934: } else {
1935: PetscArrayzero(v,PetscInt64Mult(m,n));
1936: }
1937: MatDenseRestoreArrayWrite(A,&v);
1938: return 0;
1939: }
1941: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1942: {
1943: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1944: PetscInt m = l->lda, n = A->cmap->n, i,j;
1945: PetscScalar *slot,*bb,*v;
1946: const PetscScalar *xx;
1948: if (PetscDefined(USE_DEBUG)) {
1949: for (i=0; i<N; i++) {
1952: }
1953: }
1954: if (!N) return 0;
1956: /* fix right hand side if needed */
1957: if (x && b) {
1958: VecGetArrayRead(x,&xx);
1959: VecGetArray(b,&bb);
1960: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1961: VecRestoreArrayRead(x,&xx);
1962: VecRestoreArray(b,&bb);
1963: }
1965: MatDenseGetArray(A,&v);
1966: for (i=0; i<N; i++) {
1967: slot = v + rows[i];
1968: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1969: }
1970: if (diag != 0.0) {
1972: for (i=0; i<N; i++) {
1973: slot = v + (m+1)*rows[i];
1974: *slot = diag;
1975: }
1976: }
1977: MatDenseRestoreArray(A,&v);
1978: return 0;
1979: }
1981: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1982: {
1983: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1985: *lda = mat->lda;
1986: return 0;
1987: }
1989: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
1990: {
1991: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1994: *array = mat->v;
1995: return 0;
1996: }
1998: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
1999: {
2000: if (array) *array = NULL;
2001: return 0;
2002: }
2004: /*@
2005: MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
2007: Not collective
2009: Input Parameter:
2010: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2012: Output Parameter:
2013: . lda - the leading dimension
2015: Level: intermediate
2017: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
2018: @*/
2019: PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
2020: {
2023: MatCheckPreallocated(A,1);
2024: PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
2025: return 0;
2026: }
2028: /*@
2029: MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix
2031: Not collective
2033: Input Parameters:
2034: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2035: - lda - the leading dimension
2037: Level: intermediate
2039: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
2040: @*/
2041: PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda)
2042: {
2044: PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
2045: return 0;
2046: }
2048: /*@C
2049: MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored
2051: Logically Collective on Mat
2053: Input Parameter:
2054: . mat - a dense matrix
2056: Output Parameter:
2057: . array - pointer to the data
2059: Level: intermediate
2061: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2062: @*/
2063: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
2064: {
2067: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
2068: return 0;
2069: }
2071: /*@C
2072: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
2074: Logically Collective on Mat
2076: Input Parameters:
2077: + mat - a dense matrix
2078: - array - pointer to the data
2080: Level: intermediate
2082: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2083: @*/
2084: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
2085: {
2088: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
2089: PetscObjectStateIncrease((PetscObject)A);
2090: #if defined(PETSC_HAVE_CUDA)
2091: A->offloadmask = PETSC_OFFLOAD_CPU;
2092: #endif
2093: return 0;
2094: }
2096: /*@C
2097: MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored
2099: Not Collective
2101: Input Parameter:
2102: . mat - a dense matrix
2104: Output Parameter:
2105: . array - pointer to the data
2107: Level: intermediate
2109: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2110: @*/
2111: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
2112: {
2115: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
2116: return 0;
2117: }
2119: /*@C
2120: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()
2122: Not Collective
2124: Input Parameters:
2125: + mat - a dense matrix
2126: - array - pointer to the data
2128: Level: intermediate
2130: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2131: @*/
2132: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
2133: {
2136: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
2137: return 0;
2138: }
2140: /*@C
2141: MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored
2143: Not Collective
2145: Input Parameter:
2146: . mat - a dense matrix
2148: Output Parameter:
2149: . array - pointer to the data
2151: Level: intermediate
2153: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2154: @*/
2155: PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array)
2156: {
2159: PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
2160: return 0;
2161: }
2163: /*@C
2164: MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()
2166: Not Collective
2168: Input Parameters:
2169: + mat - a dense matrix
2170: - array - pointer to the data
2172: Level: intermediate
2174: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2175: @*/
2176: PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2177: {
2180: PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2181: PetscObjectStateIncrease((PetscObject)A);
2182: #if defined(PETSC_HAVE_CUDA)
2183: A->offloadmask = PETSC_OFFLOAD_CPU;
2184: #endif
2185: return 0;
2186: }
2188: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
2189: {
2190: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2191: PetscInt i,j,nrows,ncols,ldb;
2192: const PetscInt *irow,*icol;
2193: PetscScalar *av,*bv,*v = mat->v;
2194: Mat newmat;
2196: ISGetIndices(isrow,&irow);
2197: ISGetIndices(iscol,&icol);
2198: ISGetLocalSize(isrow,&nrows);
2199: ISGetLocalSize(iscol,&ncols);
2201: /* Check submatrixcall */
2202: if (scall == MAT_REUSE_MATRIX) {
2203: PetscInt n_cols,n_rows;
2204: MatGetSize(*B,&n_rows,&n_cols);
2205: if (n_rows != nrows || n_cols != ncols) {
2206: /* resize the result matrix to match number of requested rows/columns */
2207: MatSetSizes(*B,nrows,ncols,nrows,ncols);
2208: }
2209: newmat = *B;
2210: } else {
2211: /* Create and fill new matrix */
2212: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2213: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2214: MatSetType(newmat,((PetscObject)A)->type_name);
2215: MatSeqDenseSetPreallocation(newmat,NULL);
2216: }
2218: /* Now extract the data pointers and do the copy,column at a time */
2219: MatDenseGetArray(newmat,&bv);
2220: MatDenseGetLDA(newmat,&ldb);
2221: for (i=0; i<ncols; i++) {
2222: av = v + mat->lda*icol[i];
2223: for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2224: bv += ldb;
2225: }
2226: MatDenseRestoreArray(newmat,&bv);
2228: /* Assemble the matrices so that the correct flags are set */
2229: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2230: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
2232: /* Free work space */
2233: ISRestoreIndices(isrow,&irow);
2234: ISRestoreIndices(iscol,&icol);
2235: *B = newmat;
2236: return 0;
2237: }
2239: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2240: {
2241: PetscInt i;
2243: if (scall == MAT_INITIAL_MATRIX) {
2244: PetscCalloc1(n,B);
2245: }
2247: for (i=0; i<n; i++) {
2248: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
2249: }
2250: return 0;
2251: }
2253: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2254: {
2255: return 0;
2256: }
2258: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2259: {
2260: return 0;
2261: }
2263: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2264: {
2265: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2266: const PetscScalar *va;
2267: PetscScalar *vb;
2268: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
2270: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2271: if (A->ops->copy != B->ops->copy) {
2272: MatCopy_Basic(A,B,str);
2273: return 0;
2274: }
2276: MatDenseGetArrayRead(A,&va);
2277: MatDenseGetArray(B,&vb);
2278: if (lda1>m || lda2>m) {
2279: for (j=0; j<n; j++) {
2280: PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2281: }
2282: } else {
2283: PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2284: }
2285: MatDenseRestoreArray(B,&vb);
2286: MatDenseRestoreArrayRead(A,&va);
2287: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2288: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2289: return 0;
2290: }
2292: PetscErrorCode MatSetUp_SeqDense(Mat A)
2293: {
2294: PetscLayoutSetUp(A->rmap);
2295: PetscLayoutSetUp(A->cmap);
2296: if (!A->preallocated) {
2297: MatSeqDenseSetPreallocation(A,NULL);
2298: }
2299: return 0;
2300: }
2302: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2303: {
2304: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
2305: PetscInt i,nz = A->rmap->n*A->cmap->n;
2306: PetscInt min = PetscMin(A->rmap->n,A->cmap->n);
2307: PetscScalar *aa;
2309: MatDenseGetArray(A,&aa);
2310: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2311: MatDenseRestoreArray(A,&aa);
2312: if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2313: return 0;
2314: }
2316: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2317: {
2318: PetscInt i,nz = A->rmap->n*A->cmap->n;
2319: PetscScalar *aa;
2321: MatDenseGetArray(A,&aa);
2322: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2323: MatDenseRestoreArray(A,&aa);
2324: return 0;
2325: }
2327: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2328: {
2329: PetscInt i,nz = A->rmap->n*A->cmap->n;
2330: PetscScalar *aa;
2332: MatDenseGetArray(A,&aa);
2333: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2334: MatDenseRestoreArray(A,&aa);
2335: return 0;
2336: }
2338: /* ----------------------------------------------------------------*/
2339: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2340: {
2341: PetscInt m=A->rmap->n,n=B->cmap->n;
2342: PetscBool cisdense;
2344: MatSetSizes(C,m,n,m,n);
2345: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2346: if (!cisdense) {
2347: PetscBool flg;
2349: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2350: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2351: }
2352: MatSetUp(C);
2353: return 0;
2354: }
2356: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2357: {
2358: Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2359: PetscBLASInt m,n,k;
2360: const PetscScalar *av,*bv;
2361: PetscScalar *cv;
2362: PetscScalar _DOne=1.0,_DZero=0.0;
2364: PetscBLASIntCast(C->rmap->n,&m);
2365: PetscBLASIntCast(C->cmap->n,&n);
2366: PetscBLASIntCast(A->cmap->n,&k);
2367: if (!m || !n || !k) return 0;
2368: MatDenseGetArrayRead(A,&av);
2369: MatDenseGetArrayRead(B,&bv);
2370: MatDenseGetArrayWrite(C,&cv);
2371: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2372: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2373: MatDenseRestoreArrayRead(A,&av);
2374: MatDenseRestoreArrayRead(B,&bv);
2375: MatDenseRestoreArrayWrite(C,&cv);
2376: return 0;
2377: }
2379: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2380: {
2381: PetscInt m=A->rmap->n,n=B->rmap->n;
2382: PetscBool cisdense;
2384: MatSetSizes(C,m,n,m,n);
2385: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2386: if (!cisdense) {
2387: PetscBool flg;
2389: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2390: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2391: }
2392: MatSetUp(C);
2393: return 0;
2394: }
2396: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2397: {
2398: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2399: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2400: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2401: const PetscScalar *av,*bv;
2402: PetscScalar *cv;
2403: PetscBLASInt m,n,k;
2404: PetscScalar _DOne=1.0,_DZero=0.0;
2406: PetscBLASIntCast(C->rmap->n,&m);
2407: PetscBLASIntCast(C->cmap->n,&n);
2408: PetscBLASIntCast(A->cmap->n,&k);
2409: if (!m || !n || !k) return 0;
2410: MatDenseGetArrayRead(A,&av);
2411: MatDenseGetArrayRead(B,&bv);
2412: MatDenseGetArrayWrite(C,&cv);
2413: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2414: MatDenseRestoreArrayRead(A,&av);
2415: MatDenseRestoreArrayRead(B,&bv);
2416: MatDenseRestoreArrayWrite(C,&cv);
2417: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2418: return 0;
2419: }
2421: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2422: {
2423: PetscInt m=A->cmap->n,n=B->cmap->n;
2424: PetscBool cisdense;
2426: MatSetSizes(C,m,n,m,n);
2427: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2428: if (!cisdense) {
2429: PetscBool flg;
2431: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2432: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2433: }
2434: MatSetUp(C);
2435: return 0;
2436: }
2438: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2439: {
2440: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2441: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2442: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2443: const PetscScalar *av,*bv;
2444: PetscScalar *cv;
2445: PetscBLASInt m,n,k;
2446: PetscScalar _DOne=1.0,_DZero=0.0;
2448: PetscBLASIntCast(C->rmap->n,&m);
2449: PetscBLASIntCast(C->cmap->n,&n);
2450: PetscBLASIntCast(A->rmap->n,&k);
2451: if (!m || !n || !k) return 0;
2452: MatDenseGetArrayRead(A,&av);
2453: MatDenseGetArrayRead(B,&bv);
2454: MatDenseGetArrayWrite(C,&cv);
2455: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2456: MatDenseRestoreArrayRead(A,&av);
2457: MatDenseRestoreArrayRead(B,&bv);
2458: MatDenseRestoreArrayWrite(C,&cv);
2459: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2460: return 0;
2461: }
2463: /* ----------------------------------------------- */
2464: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2465: {
2466: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2467: C->ops->productsymbolic = MatProductSymbolic_AB;
2468: return 0;
2469: }
2471: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2472: {
2473: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2474: C->ops->productsymbolic = MatProductSymbolic_AtB;
2475: return 0;
2476: }
2478: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2479: {
2480: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2481: C->ops->productsymbolic = MatProductSymbolic_ABt;
2482: return 0;
2483: }
2485: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2486: {
2487: Mat_Product *product = C->product;
2489: switch (product->type) {
2490: case MATPRODUCT_AB:
2491: MatProductSetFromOptions_SeqDense_AB(C);
2492: break;
2493: case MATPRODUCT_AtB:
2494: MatProductSetFromOptions_SeqDense_AtB(C);
2495: break;
2496: case MATPRODUCT_ABt:
2497: MatProductSetFromOptions_SeqDense_ABt(C);
2498: break;
2499: default:
2500: break;
2501: }
2502: return 0;
2503: }
2504: /* ----------------------------------------------- */
2506: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2507: {
2508: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2509: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2510: PetscScalar *x;
2511: const PetscScalar *aa;
2514: VecGetArray(v,&x);
2515: VecGetLocalSize(v,&p);
2516: MatDenseGetArrayRead(A,&aa);
2518: for (i=0; i<m; i++) {
2519: x[i] = aa[i]; if (idx) idx[i] = 0;
2520: for (j=1; j<n; j++) {
2521: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2522: }
2523: }
2524: MatDenseRestoreArrayRead(A,&aa);
2525: VecRestoreArray(v,&x);
2526: return 0;
2527: }
2529: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2530: {
2531: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2532: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2533: PetscScalar *x;
2534: PetscReal atmp;
2535: const PetscScalar *aa;
2538: VecGetArray(v,&x);
2539: VecGetLocalSize(v,&p);
2540: MatDenseGetArrayRead(A,&aa);
2542: for (i=0; i<m; i++) {
2543: x[i] = PetscAbsScalar(aa[i]);
2544: for (j=1; j<n; j++) {
2545: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2546: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2547: }
2548: }
2549: MatDenseRestoreArrayRead(A,&aa);
2550: VecRestoreArray(v,&x);
2551: return 0;
2552: }
2554: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2555: {
2556: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2557: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2558: PetscScalar *x;
2559: const PetscScalar *aa;
2562: MatDenseGetArrayRead(A,&aa);
2563: VecGetArray(v,&x);
2564: VecGetLocalSize(v,&p);
2566: for (i=0; i<m; i++) {
2567: x[i] = aa[i]; if (idx) idx[i] = 0;
2568: for (j=1; j<n; j++) {
2569: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2570: }
2571: }
2572: VecRestoreArray(v,&x);
2573: MatDenseRestoreArrayRead(A,&aa);
2574: return 0;
2575: }
2577: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2578: {
2579: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2580: PetscScalar *x;
2581: const PetscScalar *aa;
2584: MatDenseGetArrayRead(A,&aa);
2585: VecGetArray(v,&x);
2586: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2587: VecRestoreArray(v,&x);
2588: MatDenseRestoreArrayRead(A,&aa);
2589: return 0;
2590: }
2592: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions)
2593: {
2594: PetscInt i,j,m,n;
2595: const PetscScalar *a;
2597: MatGetSize(A,&m,&n);
2598: PetscArrayzero(reductions,n);
2599: MatDenseGetArrayRead(A,&a);
2600: if (type == NORM_2) {
2601: for (i=0; i<n; i++) {
2602: for (j=0; j<m; j++) {
2603: reductions[i] += PetscAbsScalar(a[j]*a[j]);
2604: }
2605: a += m;
2606: }
2607: } else if (type == NORM_1) {
2608: for (i=0; i<n; i++) {
2609: for (j=0; j<m; j++) {
2610: reductions[i] += PetscAbsScalar(a[j]);
2611: }
2612: a += m;
2613: }
2614: } else if (type == NORM_INFINITY) {
2615: for (i=0; i<n; i++) {
2616: for (j=0; j<m; j++) {
2617: reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]);
2618: }
2619: a += m;
2620: }
2621: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2622: for (i=0; i<n; i++) {
2623: for (j=0; j<m; j++) {
2624: reductions[i] += PetscRealPart(a[j]);
2625: }
2626: a += m;
2627: }
2628: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2629: for (i=0; i<n; i++) {
2630: for (j=0; j<m; j++) {
2631: reductions[i] += PetscImaginaryPart(a[j]);
2632: }
2633: a += m;
2634: }
2635: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type");
2636: MatDenseRestoreArrayRead(A,&a);
2637: if (type == NORM_2) {
2638: for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2639: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2640: for (i=0; i<n; i++) reductions[i] /= m;
2641: }
2642: return 0;
2643: }
2645: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2646: {
2647: PetscScalar *a;
2648: PetscInt lda,m,n,i,j;
2650: MatGetSize(x,&m,&n);
2651: MatDenseGetLDA(x,&lda);
2652: MatDenseGetArray(x,&a);
2653: for (j=0; j<n; j++) {
2654: for (i=0; i<m; i++) {
2655: PetscRandomGetValue(rctx,a+j*lda+i);
2656: }
2657: }
2658: MatDenseRestoreArray(x,&a);
2659: return 0;
2660: }
2662: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2663: {
2664: *missing = PETSC_FALSE;
2665: return 0;
2666: }
2668: /* vals is not const */
2669: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2670: {
2671: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2672: PetscScalar *v;
2675: MatDenseGetArray(A,&v);
2676: *vals = v+col*a->lda;
2677: MatDenseRestoreArray(A,&v);
2678: return 0;
2679: }
2681: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2682: {
2683: *vals = NULL; /* user cannot accidentally use the array later */
2684: return 0;
2685: }
2687: /* -------------------------------------------------------------------*/
2688: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2689: MatGetRow_SeqDense,
2690: MatRestoreRow_SeqDense,
2691: MatMult_SeqDense,
2692: /* 4*/ MatMultAdd_SeqDense,
2693: MatMultTranspose_SeqDense,
2694: MatMultTransposeAdd_SeqDense,
2695: NULL,
2696: NULL,
2697: NULL,
2698: /* 10*/ NULL,
2699: MatLUFactor_SeqDense,
2700: MatCholeskyFactor_SeqDense,
2701: MatSOR_SeqDense,
2702: MatTranspose_SeqDense,
2703: /* 15*/ MatGetInfo_SeqDense,
2704: MatEqual_SeqDense,
2705: MatGetDiagonal_SeqDense,
2706: MatDiagonalScale_SeqDense,
2707: MatNorm_SeqDense,
2708: /* 20*/ MatAssemblyBegin_SeqDense,
2709: MatAssemblyEnd_SeqDense,
2710: MatSetOption_SeqDense,
2711: MatZeroEntries_SeqDense,
2712: /* 24*/ MatZeroRows_SeqDense,
2713: NULL,
2714: NULL,
2715: NULL,
2716: NULL,
2717: /* 29*/ MatSetUp_SeqDense,
2718: NULL,
2719: NULL,
2720: NULL,
2721: NULL,
2722: /* 34*/ MatDuplicate_SeqDense,
2723: NULL,
2724: NULL,
2725: NULL,
2726: NULL,
2727: /* 39*/ MatAXPY_SeqDense,
2728: MatCreateSubMatrices_SeqDense,
2729: NULL,
2730: MatGetValues_SeqDense,
2731: MatCopy_SeqDense,
2732: /* 44*/ MatGetRowMax_SeqDense,
2733: MatScale_SeqDense,
2734: MatShift_Basic,
2735: NULL,
2736: MatZeroRowsColumns_SeqDense,
2737: /* 49*/ MatSetRandom_SeqDense,
2738: NULL,
2739: NULL,
2740: NULL,
2741: NULL,
2742: /* 54*/ NULL,
2743: NULL,
2744: NULL,
2745: NULL,
2746: NULL,
2747: /* 59*/ MatCreateSubMatrix_SeqDense,
2748: MatDestroy_SeqDense,
2749: MatView_SeqDense,
2750: NULL,
2751: NULL,
2752: /* 64*/ NULL,
2753: NULL,
2754: NULL,
2755: NULL,
2756: NULL,
2757: /* 69*/ MatGetRowMaxAbs_SeqDense,
2758: NULL,
2759: NULL,
2760: NULL,
2761: NULL,
2762: /* 74*/ NULL,
2763: NULL,
2764: NULL,
2765: NULL,
2766: NULL,
2767: /* 79*/ NULL,
2768: NULL,
2769: NULL,
2770: NULL,
2771: /* 83*/ MatLoad_SeqDense,
2772: MatIsSymmetric_SeqDense,
2773: MatIsHermitian_SeqDense,
2774: NULL,
2775: NULL,
2776: NULL,
2777: /* 89*/ NULL,
2778: NULL,
2779: MatMatMultNumeric_SeqDense_SeqDense,
2780: NULL,
2781: NULL,
2782: /* 94*/ NULL,
2783: NULL,
2784: NULL,
2785: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2786: NULL,
2787: /* 99*/ MatProductSetFromOptions_SeqDense,
2788: NULL,
2789: NULL,
2790: MatConjugate_SeqDense,
2791: NULL,
2792: /*104*/ NULL,
2793: MatRealPart_SeqDense,
2794: MatImaginaryPart_SeqDense,
2795: NULL,
2796: NULL,
2797: /*109*/ NULL,
2798: NULL,
2799: MatGetRowMin_SeqDense,
2800: MatGetColumnVector_SeqDense,
2801: MatMissingDiagonal_SeqDense,
2802: /*114*/ NULL,
2803: NULL,
2804: NULL,
2805: NULL,
2806: NULL,
2807: /*119*/ NULL,
2808: NULL,
2809: NULL,
2810: NULL,
2811: NULL,
2812: /*124*/ NULL,
2813: MatGetColumnReductions_SeqDense,
2814: NULL,
2815: NULL,
2816: NULL,
2817: /*129*/ NULL,
2818: NULL,
2819: NULL,
2820: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2821: NULL,
2822: /*134*/ NULL,
2823: NULL,
2824: NULL,
2825: NULL,
2826: NULL,
2827: /*139*/ NULL,
2828: NULL,
2829: NULL,
2830: NULL,
2831: NULL,
2832: MatCreateMPIMatConcatenateSeqMat_SeqDense,
2833: /*145*/ NULL,
2834: NULL,
2835: NULL
2836: };
2838: /*@C
2839: MatCreateSeqDense - Creates a sequential dense matrix that
2840: is stored in column major order (the usual Fortran 77 manner). Many
2841: of the matrix operations use the BLAS and LAPACK routines.
2843: Collective
2845: Input Parameters:
2846: + comm - MPI communicator, set to PETSC_COMM_SELF
2847: . m - number of rows
2848: . n - number of columns
2849: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2850: to control all matrix memory allocation.
2852: Output Parameter:
2853: . A - the matrix
2855: Notes:
2856: The data input variable is intended primarily for Fortran programmers
2857: who wish to allocate their own matrix memory space. Most users should
2858: set data=NULL.
2860: Level: intermediate
2862: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2863: @*/
2864: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2865: {
2866: MatCreate(comm,A);
2867: MatSetSizes(*A,m,n,m,n);
2868: MatSetType(*A,MATSEQDENSE);
2869: MatSeqDenseSetPreallocation(*A,data);
2870: return 0;
2871: }
2873: /*@C
2874: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2876: Collective
2878: Input Parameters:
2879: + B - the matrix
2880: - data - the array (or NULL)
2882: Notes:
2883: The data input variable is intended primarily for Fortran programmers
2884: who wish to allocate their own matrix memory space. Most users should
2885: need not call this routine.
2887: Level: intermediate
2889: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
2891: @*/
2892: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2893: {
2895: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2896: return 0;
2897: }
2899: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2900: {
2901: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2904: B->preallocated = PETSC_TRUE;
2906: PetscLayoutSetUp(B->rmap);
2907: PetscLayoutSetUp(B->cmap);
2909: if (b->lda <= 0) b->lda = B->rmap->n;
2911: if (!data) { /* petsc-allocated storage */
2912: if (!b->user_alloc) PetscFree(b->v);
2913: PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
2914: PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));
2916: b->user_alloc = PETSC_FALSE;
2917: } else { /* user-allocated storage */
2918: if (!b->user_alloc) PetscFree(b->v);
2919: b->v = data;
2920: b->user_alloc = PETSC_TRUE;
2921: }
2922: B->assembled = PETSC_TRUE;
2923: return 0;
2924: }
2926: #if defined(PETSC_HAVE_ELEMENTAL)
2927: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2928: {
2929: Mat mat_elemental;
2930: const PetscScalar *array;
2931: PetscScalar *v_colwise;
2932: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2934: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2935: MatDenseGetArrayRead(A,&array);
2936: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2937: k = 0;
2938: for (j=0; j<N; j++) {
2939: cols[j] = j;
2940: for (i=0; i<M; i++) {
2941: v_colwise[j*M+i] = array[k++];
2942: }
2943: }
2944: for (i=0; i<M; i++) {
2945: rows[i] = i;
2946: }
2947: MatDenseRestoreArrayRead(A,&array);
2949: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2950: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2951: MatSetType(mat_elemental,MATELEMENTAL);
2952: MatSetUp(mat_elemental);
2954: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2955: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2956: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2957: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2958: PetscFree3(v_colwise,rows,cols);
2960: if (reuse == MAT_INPLACE_MATRIX) {
2961: MatHeaderReplace(A,&mat_elemental);
2962: } else {
2963: *newmat = mat_elemental;
2964: }
2965: return 0;
2966: }
2967: #endif
2969: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2970: {
2971: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2972: PetscBool data;
2974: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2977: b->lda = lda;
2978: return 0;
2979: }
2981: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2982: {
2983: PetscMPIInt size;
2985: MPI_Comm_size(comm,&size);
2986: if (size == 1) {
2987: if (scall == MAT_INITIAL_MATRIX) {
2988: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2989: } else {
2990: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2991: }
2992: } else {
2993: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2994: }
2995: return 0;
2996: }
2998: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2999: {
3000: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3004: if (!a->cvec) {
3005: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3006: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3007: }
3008: a->vecinuse = col + 1;
3009: MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3010: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3011: *v = a->cvec;
3012: return 0;
3013: }
3015: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3016: {
3017: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3021: a->vecinuse = 0;
3022: MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3023: VecResetArray(a->cvec);
3024: if (v) *v = NULL;
3025: return 0;
3026: }
3028: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3029: {
3030: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3034: if (!a->cvec) {
3035: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3036: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3037: }
3038: a->vecinuse = col + 1;
3039: MatDenseGetArrayRead(A,&a->ptrinuse);
3040: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3041: VecLockReadPush(a->cvec);
3042: *v = a->cvec;
3043: return 0;
3044: }
3046: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3047: {
3048: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3052: a->vecinuse = 0;
3053: MatDenseRestoreArrayRead(A,&a->ptrinuse);
3054: VecLockReadPop(a->cvec);
3055: VecResetArray(a->cvec);
3056: if (v) *v = NULL;
3057: return 0;
3058: }
3060: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3061: {
3062: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3066: if (!a->cvec) {
3067: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3068: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3069: }
3070: a->vecinuse = col + 1;
3071: MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3072: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3073: *v = a->cvec;
3074: return 0;
3075: }
3077: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3078: {
3079: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3083: a->vecinuse = 0;
3084: MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3085: VecResetArray(a->cvec);
3086: if (v) *v = NULL;
3087: return 0;
3088: }
3090: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3091: {
3092: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3096: if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3097: MatDestroy(&a->cmat);
3098: }
3099: if (!a->cmat) {
3100: MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,a->v+(size_t)cbegin*a->lda,&a->cmat);
3101: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3102: } else {
3103: MatDensePlaceArray(a->cmat,a->v+(size_t)cbegin*a->lda);
3104: }
3105: MatDenseSetLDA(a->cmat,a->lda);
3106: a->matinuse = cbegin + 1;
3107: *v = a->cmat;
3108: #if defined(PETSC_HAVE_CUDA)
3109: A->offloadmask = PETSC_OFFLOAD_CPU;
3110: #endif
3111: return 0;
3112: }
3114: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3115: {
3116: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3121: a->matinuse = 0;
3122: MatDenseResetArray(a->cmat);
3123: *v = NULL;
3124: return 0;
3125: }
3127: /*MC
3128: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3130: Options Database Keys:
3131: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3133: Level: beginner
3135: .seealso: MatCreateSeqDense()
3137: M*/
3138: PetscErrorCode MatCreate_SeqDense(Mat B)
3139: {
3140: Mat_SeqDense *b;
3141: PetscMPIInt size;
3143: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3146: PetscNewLog(B,&b);
3147: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3148: B->data = (void*)b;
3150: b->roworiented = PETSC_TRUE;
3152: PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3153: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3154: PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3155: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3156: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3157: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3158: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3159: PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3160: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3161: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3162: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3163: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3164: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3165: #if defined(PETSC_HAVE_ELEMENTAL)
3166: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3167: #endif
3168: #if defined(PETSC_HAVE_SCALAPACK)
3169: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3170: #endif
3171: #if defined(PETSC_HAVE_CUDA)
3172: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3173: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3174: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3175: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3176: #endif
3177: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3178: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3179: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3180: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3181: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3183: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3184: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3185: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3186: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3187: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3188: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3189: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3190: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3191: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3192: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3193: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3194: return 0;
3195: }
3197: /*@C
3198: MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding.
3200: Not Collective
3202: Input Parameters:
3203: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3204: - col - column index
3206: Output Parameter:
3207: . vals - pointer to the data
3209: Level: intermediate
3211: .seealso: MatDenseRestoreColumn()
3212: @*/
3213: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3214: {
3218: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3219: return 0;
3220: }
3222: /*@C
3223: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3225: Not Collective
3227: Input Parameter:
3228: . mat - a MATSEQDENSE or MATMPIDENSE matrix
3230: Output Parameter:
3231: . vals - pointer to the data
3233: Level: intermediate
3235: .seealso: MatDenseGetColumn()
3236: @*/
3237: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3238: {
3241: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3242: return 0;
3243: }
3245: /*@
3246: MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3248: Collective
3250: Input Parameters:
3251: + mat - the Mat object
3252: - col - the column index
3254: Output Parameter:
3255: . v - the vector
3257: Notes:
3258: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3259: Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.
3261: Level: intermediate
3263: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3264: @*/
3265: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3266: {
3273: PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3274: return 0;
3275: }
3277: /*@
3278: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3280: Collective
3282: Input Parameters:
3283: + mat - the Mat object
3284: . col - the column index
3285: - v - the Vec object
3287: Level: intermediate
3289: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3290: @*/
3291: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3292: {
3298: PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3299: return 0;
3300: }
3302: /*@
3303: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3305: Collective
3307: Input Parameters:
3308: + mat - the Mat object
3309: - col - the column index
3311: Output Parameter:
3312: . v - the vector
3314: Notes:
3315: The vector is owned by PETSc and users cannot modify it.
3316: Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3317: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.
3319: Level: intermediate
3321: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3322: @*/
3323: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3324: {
3331: PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3332: return 0;
3333: }
3335: /*@
3336: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3338: Collective
3340: Input Parameters:
3341: + mat - the Mat object
3342: . col - the column index
3343: - v - the Vec object
3345: Level: intermediate
3347: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3348: @*/
3349: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3350: {
3356: PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3357: return 0;
3358: }
3360: /*@
3361: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3363: Collective
3365: Input Parameters:
3366: + mat - the Mat object
3367: - col - the column index
3369: Output Parameter:
3370: . v - the vector
3372: Notes:
3373: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3374: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.
3376: Level: intermediate
3378: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3379: @*/
3380: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3381: {
3388: PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3389: return 0;
3390: }
3392: /*@
3393: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3395: Collective
3397: Input Parameters:
3398: + mat - the Mat object
3399: . col - the column index
3400: - v - the Vec object
3402: Level: intermediate
3404: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3405: @*/
3406: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3407: {
3413: PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3414: return 0;
3415: }
3417: /*@
3418: MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3420: Collective
3422: Input Parameters:
3423: + mat - the Mat object
3424: . cbegin - the first index in the block
3425: - cend - the last index in the block
3427: Output Parameter:
3428: . v - the matrix
3430: Notes:
3431: The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3433: Level: intermediate
3435: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3436: @*/
3437: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3438: {
3447: PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3448: return 0;
3449: }
3451: /*@
3452: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3454: Collective
3456: Input Parameters:
3457: + mat - the Mat object
3458: - v - the Mat object
3460: Level: intermediate
3462: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3463: @*/
3464: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3465: {
3469: PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3470: return 0;
3471: }