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: /* ---------------------------------------------------------------*/
763: /* COMMENT: I have chosen to hide row permutation in the pivots,
764: rather than put it in the Mat->row slot.*/
765: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
766: {
767: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
768: PetscBLASInt n,m,info;
770: PetscBLASIntCast(A->cmap->n,&n);
771: PetscBLASIntCast(A->rmap->n,&m);
772: if (!mat->pivots) {
773: PetscMalloc1(A->rmap->n,&mat->pivots);
774: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
775: }
776: if (!A->rmap->n || !A->cmap->n) return 0;
777: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
778: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
779: PetscFPTrapPop();
784: A->ops->solve = MatSolve_SeqDense_LU;
785: A->ops->matsolve = MatMatSolve_SeqDense_LU;
786: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
787: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
788: A->factortype = MAT_FACTOR_LU;
790: PetscFree(A->solvertype);
791: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
793: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
794: return 0;
795: }
797: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
798: {
799: MatFactorInfo info;
801: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
802: (*fact->ops->lufactor)(fact,NULL,NULL,&info);
803: return 0;
804: }
806: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
807: {
808: fact->preallocated = PETSC_TRUE;
809: fact->assembled = PETSC_TRUE;
810: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
811: return 0;
812: }
814: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
815: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
816: {
817: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
818: PetscBLASInt info,n;
820: PetscBLASIntCast(A->cmap->n,&n);
821: if (!A->rmap->n || !A->cmap->n) return 0;
822: if (A->spd) {
823: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
824: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
825: PetscFPTrapPop();
826: #if defined(PETSC_USE_COMPLEX)
827: } else if (A->hermitian) {
828: if (!mat->pivots) {
829: PetscMalloc1(A->rmap->n,&mat->pivots);
830: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
831: }
832: if (!mat->fwork) {
833: PetscScalar dummy;
835: mat->lfwork = -1;
836: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
837: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
838: PetscFPTrapPop();
839: mat->lfwork = (PetscInt)PetscRealPart(dummy);
840: PetscMalloc1(mat->lfwork,&mat->fwork);
841: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
842: }
843: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
844: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
845: PetscFPTrapPop();
846: #endif
847: } else { /* symmetric case */
848: if (!mat->pivots) {
849: PetscMalloc1(A->rmap->n,&mat->pivots);
850: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
851: }
852: if (!mat->fwork) {
853: PetscScalar dummy;
855: mat->lfwork = -1;
856: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
857: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
858: PetscFPTrapPop();
859: mat->lfwork = (PetscInt)PetscRealPart(dummy);
860: PetscMalloc1(mat->lfwork,&mat->fwork);
861: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
862: }
863: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
864: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
865: PetscFPTrapPop();
866: }
869: A->ops->solve = MatSolve_SeqDense_Cholesky;
870: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
871: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
872: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
873: A->factortype = MAT_FACTOR_CHOLESKY;
875: PetscFree(A->solvertype);
876: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
878: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
879: return 0;
880: }
882: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
883: {
884: MatFactorInfo info;
886: info.fill = 1.0;
888: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
889: (*fact->ops->choleskyfactor)(fact,NULL,&info);
890: return 0;
891: }
893: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
894: {
895: fact->assembled = PETSC_TRUE;
896: fact->preallocated = PETSC_TRUE;
897: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
898: return 0;
899: }
901: PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
902: {
903: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
904: PetscBLASInt n,m,info, min, max;
906: PetscBLASIntCast(A->cmap->n,&n);
907: PetscBLASIntCast(A->rmap->n,&m);
908: max = PetscMax(m, n);
909: min = PetscMin(m, n);
910: if (!mat->tau) {
911: PetscMalloc1(min,&mat->tau);
912: PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
913: }
914: if (!mat->pivots) {
915: PetscMalloc1(n,&mat->pivots);
916: PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar));
917: }
918: if (!mat->qrrhs) {
919: MatCreateVecs(A, NULL, &(mat->qrrhs));
920: }
921: if (!A->rmap->n || !A->cmap->n) return 0;
922: if (!mat->fwork) {
923: PetscScalar dummy;
925: mat->lfwork = -1;
926: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
927: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
928: PetscFPTrapPop();
929: mat->lfwork = (PetscInt)PetscRealPart(dummy);
930: PetscMalloc1(mat->lfwork,&mat->fwork);
931: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
932: }
933: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
934: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
935: PetscFPTrapPop();
937: // 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
938: mat->rank = min;
940: A->ops->solve = MatSolve_SeqDense_QR;
941: A->ops->matsolve = MatMatSolve_SeqDense_QR;
942: A->factortype = MAT_FACTOR_QR;
943: if (m == n) {
944: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
945: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
946: }
948: PetscFree(A->solvertype);
949: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
951: PetscLogFlops(2.0*min*min*(max-min/3.0));
952: return 0;
953: }
955: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
956: {
957: MatFactorInfo info;
959: info.fill = 1.0;
961: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
962: PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info));
963: return 0;
964: }
966: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
967: {
968: fact->assembled = PETSC_TRUE;
969: fact->preallocated = PETSC_TRUE;
970: PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
971: return 0;
972: }
974: /* uses LAPACK */
975: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
976: {
977: MatCreate(PetscObjectComm((PetscObject)A),fact);
978: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
979: MatSetType(*fact,MATDENSE);
980: (*fact)->trivialsymbolic = PETSC_TRUE;
981: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
982: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
983: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
984: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
985: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
986: } else if (ftype == MAT_FACTOR_QR) {
987: PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
988: }
989: (*fact)->factortype = ftype;
991: PetscFree((*fact)->solvertype);
992: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
993: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU]);
994: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU]);
995: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]);
996: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC]);
997: return 0;
998: }
1000: /* ------------------------------------------------------------------*/
1001: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1002: {
1003: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1004: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
1005: const PetscScalar *b;
1006: PetscInt m = A->rmap->n,i;
1007: PetscBLASInt o = 1,bm = 0;
1009: #if defined(PETSC_HAVE_CUDA)
1011: #endif
1012: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1013: PetscBLASIntCast(m,&bm);
1014: if (flag & SOR_ZERO_INITIAL_GUESS) {
1015: /* this is a hack fix, should have another version without the second BLASdotu */
1016: VecSet(xx,zero);
1017: }
1018: VecGetArray(xx,&x);
1019: VecGetArrayRead(bb,&b);
1020: its = its*lits;
1022: while (its--) {
1023: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1024: for (i=0; i<m; i++) {
1025: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1026: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1027: }
1028: }
1029: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1030: for (i=m-1; i>=0; i--) {
1031: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1032: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1033: }
1034: }
1035: }
1036: VecRestoreArrayRead(bb,&b);
1037: VecRestoreArray(xx,&x);
1038: return 0;
1039: }
1041: /* -----------------------------------------------------------------*/
1042: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1043: {
1044: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1045: const PetscScalar *v = mat->v,*x;
1046: PetscScalar *y;
1047: PetscBLASInt m, n,_One=1;
1048: PetscScalar _DOne=1.0,_DZero=0.0;
1050: PetscBLASIntCast(A->rmap->n,&m);
1051: PetscBLASIntCast(A->cmap->n,&n);
1052: VecGetArrayRead(xx,&x);
1053: VecGetArrayWrite(yy,&y);
1054: if (!A->rmap->n || !A->cmap->n) {
1055: PetscBLASInt i;
1056: for (i=0; i<n; i++) y[i] = 0.0;
1057: } else {
1058: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1059: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1060: }
1061: VecRestoreArrayRead(xx,&x);
1062: VecRestoreArrayWrite(yy,&y);
1063: return 0;
1064: }
1066: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1067: {
1068: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1069: PetscScalar *y,_DOne=1.0,_DZero=0.0;
1070: PetscBLASInt m, n, _One=1;
1071: const PetscScalar *v = mat->v,*x;
1073: PetscBLASIntCast(A->rmap->n,&m);
1074: PetscBLASIntCast(A->cmap->n,&n);
1075: VecGetArrayRead(xx,&x);
1076: VecGetArrayWrite(yy,&y);
1077: if (!A->rmap->n || !A->cmap->n) {
1078: PetscBLASInt i;
1079: for (i=0; i<m; i++) y[i] = 0.0;
1080: } else {
1081: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1082: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1083: }
1084: VecRestoreArrayRead(xx,&x);
1085: VecRestoreArrayWrite(yy,&y);
1086: return 0;
1087: }
1089: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1090: {
1091: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1092: const PetscScalar *v = mat->v,*x;
1093: PetscScalar *y,_DOne=1.0;
1094: PetscBLASInt m, n, _One=1;
1096: PetscBLASIntCast(A->rmap->n,&m);
1097: PetscBLASIntCast(A->cmap->n,&n);
1098: VecCopy(zz,yy);
1099: if (!A->rmap->n || !A->cmap->n) return 0;
1100: VecGetArrayRead(xx,&x);
1101: VecGetArray(yy,&y);
1102: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1103: VecRestoreArrayRead(xx,&x);
1104: VecRestoreArray(yy,&y);
1105: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1106: return 0;
1107: }
1109: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1110: {
1111: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1112: const PetscScalar *v = mat->v,*x;
1113: PetscScalar *y;
1114: PetscBLASInt m, n, _One=1;
1115: PetscScalar _DOne=1.0;
1117: PetscBLASIntCast(A->rmap->n,&m);
1118: PetscBLASIntCast(A->cmap->n,&n);
1119: VecCopy(zz,yy);
1120: if (!A->rmap->n || !A->cmap->n) return 0;
1121: VecGetArrayRead(xx,&x);
1122: VecGetArray(yy,&y);
1123: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1124: VecRestoreArrayRead(xx,&x);
1125: VecRestoreArray(yy,&y);
1126: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1127: return 0;
1128: }
1130: /* -----------------------------------------------------------------*/
1131: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1132: {
1133: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1134: PetscInt i;
1136: *ncols = A->cmap->n;
1137: if (cols) {
1138: PetscMalloc1(A->cmap->n,cols);
1139: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1140: }
1141: if (vals) {
1142: const PetscScalar *v;
1144: MatDenseGetArrayRead(A,&v);
1145: PetscMalloc1(A->cmap->n,vals);
1146: v += row;
1147: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1148: MatDenseRestoreArrayRead(A,&v);
1149: }
1150: return 0;
1151: }
1153: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1154: {
1155: if (ncols) *ncols = 0;
1156: if (cols) PetscFree(*cols);
1157: if (vals) PetscFree(*vals);
1158: return 0;
1159: }
1160: /* ----------------------------------------------------------------*/
1161: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1162: {
1163: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1164: PetscScalar *av;
1165: PetscInt i,j,idx=0;
1166: #if defined(PETSC_HAVE_CUDA)
1167: PetscOffloadMask oldf;
1168: #endif
1170: MatDenseGetArray(A,&av);
1171: if (!mat->roworiented) {
1172: if (addv == INSERT_VALUES) {
1173: for (j=0; j<n; j++) {
1174: if (indexn[j] < 0) {idx += m; continue;}
1176: for (i=0; i<m; i++) {
1177: if (indexm[i] < 0) {idx++; continue;}
1179: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1180: }
1181: }
1182: } else {
1183: for (j=0; j<n; j++) {
1184: if (indexn[j] < 0) {idx += m; continue;}
1186: for (i=0; i<m; i++) {
1187: if (indexm[i] < 0) {idx++; continue;}
1189: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1190: }
1191: }
1192: }
1193: } else {
1194: if (addv == INSERT_VALUES) {
1195: for (i=0; i<m; i++) {
1196: if (indexm[i] < 0) { idx += n; continue;}
1198: for (j=0; j<n; j++) {
1199: if (indexn[j] < 0) { idx++; continue;}
1201: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1202: }
1203: }
1204: } else {
1205: for (i=0; i<m; i++) {
1206: if (indexm[i] < 0) { idx += n; continue;}
1208: for (j=0; j<n; j++) {
1209: if (indexn[j] < 0) { idx++; continue;}
1211: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1212: }
1213: }
1214: }
1215: }
1216: /* hack to prevent unneeded copy to the GPU while returning the array */
1217: #if defined(PETSC_HAVE_CUDA)
1218: oldf = A->offloadmask;
1219: A->offloadmask = PETSC_OFFLOAD_GPU;
1220: #endif
1221: MatDenseRestoreArray(A,&av);
1222: #if defined(PETSC_HAVE_CUDA)
1223: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1224: #endif
1225: return 0;
1226: }
1228: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1229: {
1230: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1231: const PetscScalar *vv;
1232: PetscInt i,j;
1234: MatDenseGetArrayRead(A,&vv);
1235: /* row-oriented output */
1236: for (i=0; i<m; i++) {
1237: if (indexm[i] < 0) {v += n;continue;}
1239: for (j=0; j<n; j++) {
1240: if (indexn[j] < 0) {v++; continue;}
1242: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1243: }
1244: }
1245: MatDenseRestoreArrayRead(A,&vv);
1246: return 0;
1247: }
1249: /* -----------------------------------------------------------------*/
1251: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1252: {
1253: PetscBool skipHeader;
1254: PetscViewerFormat format;
1255: PetscInt header[4],M,N,m,lda,i,j,k;
1256: const PetscScalar *v;
1257: PetscScalar *vwork;
1259: PetscViewerSetUp(viewer);
1260: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1261: PetscViewerGetFormat(viewer,&format);
1262: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1264: MatGetSize(mat,&M,&N);
1266: /* write matrix header */
1267: header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1268: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1269: if (!skipHeader) PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);
1271: MatGetLocalSize(mat,&m,NULL);
1272: if (format != PETSC_VIEWER_NATIVE) {
1273: PetscInt nnz = m*N, *iwork;
1274: /* store row lengths for each row */
1275: PetscMalloc1(nnz,&iwork);
1276: for (i=0; i<m; i++) iwork[i] = N;
1277: PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1278: /* store column indices (zero start index) */
1279: for (k=0, i=0; i<m; i++)
1280: for (j=0; j<N; j++, k++)
1281: iwork[k] = j;
1282: PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1283: PetscFree(iwork);
1284: }
1285: /* store matrix values as a dense matrix in row major order */
1286: PetscMalloc1(m*N,&vwork);
1287: MatDenseGetArrayRead(mat,&v);
1288: MatDenseGetLDA(mat,&lda);
1289: for (k=0, i=0; i<m; i++)
1290: for (j=0; j<N; j++, k++)
1291: vwork[k] = v[i+lda*j];
1292: MatDenseRestoreArrayRead(mat,&v);
1293: PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1294: PetscFree(vwork);
1295: return 0;
1296: }
1298: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1299: {
1300: PetscBool skipHeader;
1301: PetscInt header[4],M,N,m,nz,lda,i,j,k;
1302: PetscInt rows,cols;
1303: PetscScalar *v,*vwork;
1305: PetscViewerSetUp(viewer);
1306: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1308: if (!skipHeader) {
1309: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1311: M = header[1]; N = header[2];
1314: nz = header[3];
1316: } else {
1317: MatGetSize(mat,&M,&N);
1319: nz = MATRIX_BINARY_FORMAT_DENSE;
1320: }
1322: /* setup global sizes if not set */
1323: if (mat->rmap->N < 0) mat->rmap->N = M;
1324: if (mat->cmap->N < 0) mat->cmap->N = N;
1325: MatSetUp(mat);
1326: /* check if global sizes are correct */
1327: MatGetSize(mat,&rows,&cols);
1330: MatGetSize(mat,NULL,&N);
1331: MatGetLocalSize(mat,&m,NULL);
1332: MatDenseGetArray(mat,&v);
1333: MatDenseGetLDA(mat,&lda);
1334: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1335: PetscInt nnz = m*N;
1336: /* read in matrix values */
1337: PetscMalloc1(nnz,&vwork);
1338: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1339: /* store values in column major order */
1340: for (j=0; j<N; j++)
1341: for (i=0; i<m; i++)
1342: v[i+lda*j] = vwork[i*N+j];
1343: PetscFree(vwork);
1344: } else { /* matrix in file is sparse format */
1345: PetscInt nnz = 0, *rlens, *icols;
1346: /* read in row lengths */
1347: PetscMalloc1(m,&rlens);
1348: PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1349: for (i=0; i<m; i++) nnz += rlens[i];
1350: /* read in column indices and values */
1351: PetscMalloc2(nnz,&icols,nnz,&vwork);
1352: PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1353: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1354: /* store values in column major order */
1355: for (k=0, i=0; i<m; i++)
1356: for (j=0; j<rlens[i]; j++, k++)
1357: v[i+lda*icols[k]] = vwork[k];
1358: PetscFree(rlens);
1359: PetscFree2(icols,vwork);
1360: }
1361: MatDenseRestoreArray(mat,&v);
1362: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1363: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1364: return 0;
1365: }
1367: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1368: {
1369: PetscBool isbinary, ishdf5;
1373: /* force binary viewer to load .info file if it has not yet done so */
1374: PetscViewerSetUp(viewer);
1375: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1376: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1377: if (isbinary) {
1378: MatLoad_Dense_Binary(newMat,viewer);
1379: } else if (ishdf5) {
1380: #if defined(PETSC_HAVE_HDF5)
1381: MatLoad_Dense_HDF5(newMat,viewer);
1382: #else
1383: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1384: #endif
1385: } else {
1386: 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);
1387: }
1388: return 0;
1389: }
1391: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1392: {
1393: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1394: PetscInt i,j;
1395: const char *name;
1396: PetscScalar *v,*av;
1397: PetscViewerFormat format;
1398: #if defined(PETSC_USE_COMPLEX)
1399: PetscBool allreal = PETSC_TRUE;
1400: #endif
1402: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1403: PetscViewerGetFormat(viewer,&format);
1404: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1405: return 0; /* do nothing for now */
1406: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1407: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1408: for (i=0; i<A->rmap->n; i++) {
1409: v = av + i;
1410: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
1411: for (j=0; j<A->cmap->n; j++) {
1412: #if defined(PETSC_USE_COMPLEX)
1413: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1414: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1415: } else if (PetscRealPart(*v)) {
1416: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)PetscRealPart(*v));
1417: }
1418: #else
1419: if (*v) {
1420: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",j,(double)*v);
1421: }
1422: #endif
1423: v += a->lda;
1424: }
1425: PetscViewerASCIIPrintf(viewer,"\n");
1426: }
1427: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1428: } else {
1429: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1430: #if defined(PETSC_USE_COMPLEX)
1431: /* determine if matrix has all real values */
1432: for (j=0; j<A->cmap->n; j++) {
1433: v = av + j*a->lda;
1434: for (i=0; i<A->rmap->n; i++) {
1435: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1436: }
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 (may be NULL)
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 (may be NULL)
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 (may be NULL)
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,j;
2306: PetscInt min = PetscMin(A->rmap->n,A->cmap->n);
2307: PetscScalar *aa;
2309: MatDenseGetArray(A,&aa);
2310: for (j=0; j<A->cmap->n; j++) {
2311: for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscConj(aa[i+j*mat->lda]);
2312: }
2313: MatDenseRestoreArray(A,&aa);
2314: if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2315: return 0;
2316: }
2318: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2319: {
2320: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
2321: PetscInt i,j;
2322: PetscScalar *aa;
2324: MatDenseGetArray(A,&aa);
2325: for (j=0; j<A->cmap->n; j++) {
2326: for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscRealPart(aa[i+j*mat->lda]);
2327: }
2328: MatDenseRestoreArray(A,&aa);
2329: return 0;
2330: }
2332: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2333: {
2334: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
2335: PetscInt i,j;
2336: PetscScalar *aa;
2338: MatDenseGetArray(A,&aa);
2339: for (j=0; j<A->cmap->n; j++) {
2340: for (i=0; i<A->rmap->n; i++) aa[i+j*mat->lda] = PetscImaginaryPart(aa[i+j*mat->lda]);
2341: }
2342: MatDenseRestoreArray(A,&aa);
2343: return 0;
2344: }
2346: /* ----------------------------------------------------------------*/
2347: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2348: {
2349: PetscInt m=A->rmap->n,n=B->cmap->n;
2350: PetscBool cisdense;
2352: MatSetSizes(C,m,n,m,n);
2353: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2354: if (!cisdense) {
2355: PetscBool flg;
2357: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2358: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2359: }
2360: MatSetUp(C);
2361: return 0;
2362: }
2364: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2365: {
2366: Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2367: PetscBLASInt m,n,k;
2368: const PetscScalar *av,*bv;
2369: PetscScalar *cv;
2370: PetscScalar _DOne=1.0,_DZero=0.0;
2372: PetscBLASIntCast(C->rmap->n,&m);
2373: PetscBLASIntCast(C->cmap->n,&n);
2374: PetscBLASIntCast(A->cmap->n,&k);
2375: if (!m || !n || !k) return 0;
2376: MatDenseGetArrayRead(A,&av);
2377: MatDenseGetArrayRead(B,&bv);
2378: MatDenseGetArrayWrite(C,&cv);
2379: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2380: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2381: MatDenseRestoreArrayRead(A,&av);
2382: MatDenseRestoreArrayRead(B,&bv);
2383: MatDenseRestoreArrayWrite(C,&cv);
2384: return 0;
2385: }
2387: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2388: {
2389: PetscInt m=A->rmap->n,n=B->rmap->n;
2390: PetscBool cisdense;
2392: MatSetSizes(C,m,n,m,n);
2393: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2394: if (!cisdense) {
2395: PetscBool flg;
2397: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2398: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2399: }
2400: MatSetUp(C);
2401: return 0;
2402: }
2404: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2405: {
2406: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2407: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2408: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2409: const PetscScalar *av,*bv;
2410: PetscScalar *cv;
2411: PetscBLASInt m,n,k;
2412: PetscScalar _DOne=1.0,_DZero=0.0;
2414: PetscBLASIntCast(C->rmap->n,&m);
2415: PetscBLASIntCast(C->cmap->n,&n);
2416: PetscBLASIntCast(A->cmap->n,&k);
2417: if (!m || !n || !k) return 0;
2418: MatDenseGetArrayRead(A,&av);
2419: MatDenseGetArrayRead(B,&bv);
2420: MatDenseGetArrayWrite(C,&cv);
2421: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2422: MatDenseRestoreArrayRead(A,&av);
2423: MatDenseRestoreArrayRead(B,&bv);
2424: MatDenseRestoreArrayWrite(C,&cv);
2425: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2426: return 0;
2427: }
2429: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2430: {
2431: PetscInt m=A->cmap->n,n=B->cmap->n;
2432: PetscBool cisdense;
2434: MatSetSizes(C,m,n,m,n);
2435: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2436: if (!cisdense) {
2437: PetscBool flg;
2439: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2440: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2441: }
2442: MatSetUp(C);
2443: return 0;
2444: }
2446: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2447: {
2448: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2449: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2450: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2451: const PetscScalar *av,*bv;
2452: PetscScalar *cv;
2453: PetscBLASInt m,n,k;
2454: PetscScalar _DOne=1.0,_DZero=0.0;
2456: PetscBLASIntCast(C->rmap->n,&m);
2457: PetscBLASIntCast(C->cmap->n,&n);
2458: PetscBLASIntCast(A->rmap->n,&k);
2459: if (!m || !n || !k) return 0;
2460: MatDenseGetArrayRead(A,&av);
2461: MatDenseGetArrayRead(B,&bv);
2462: MatDenseGetArrayWrite(C,&cv);
2463: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2464: MatDenseRestoreArrayRead(A,&av);
2465: MatDenseRestoreArrayRead(B,&bv);
2466: MatDenseRestoreArrayWrite(C,&cv);
2467: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2468: return 0;
2469: }
2471: /* ----------------------------------------------- */
2472: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2473: {
2474: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2475: C->ops->productsymbolic = MatProductSymbolic_AB;
2476: return 0;
2477: }
2479: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2480: {
2481: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2482: C->ops->productsymbolic = MatProductSymbolic_AtB;
2483: return 0;
2484: }
2486: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2487: {
2488: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2489: C->ops->productsymbolic = MatProductSymbolic_ABt;
2490: return 0;
2491: }
2493: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2494: {
2495: Mat_Product *product = C->product;
2497: switch (product->type) {
2498: case MATPRODUCT_AB:
2499: MatProductSetFromOptions_SeqDense_AB(C);
2500: break;
2501: case MATPRODUCT_AtB:
2502: MatProductSetFromOptions_SeqDense_AtB(C);
2503: break;
2504: case MATPRODUCT_ABt:
2505: MatProductSetFromOptions_SeqDense_ABt(C);
2506: break;
2507: default:
2508: break;
2509: }
2510: return 0;
2511: }
2512: /* ----------------------------------------------- */
2514: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2515: {
2516: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2517: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2518: PetscScalar *x;
2519: const PetscScalar *aa;
2522: VecGetArray(v,&x);
2523: VecGetLocalSize(v,&p);
2524: MatDenseGetArrayRead(A,&aa);
2526: for (i=0; i<m; i++) {
2527: x[i] = aa[i]; if (idx) idx[i] = 0;
2528: for (j=1; j<n; j++) {
2529: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2530: }
2531: }
2532: MatDenseRestoreArrayRead(A,&aa);
2533: VecRestoreArray(v,&x);
2534: return 0;
2535: }
2537: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2538: {
2539: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2540: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2541: PetscScalar *x;
2542: PetscReal atmp;
2543: const PetscScalar *aa;
2546: VecGetArray(v,&x);
2547: VecGetLocalSize(v,&p);
2548: MatDenseGetArrayRead(A,&aa);
2550: for (i=0; i<m; i++) {
2551: x[i] = PetscAbsScalar(aa[i]);
2552: for (j=1; j<n; j++) {
2553: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2554: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2555: }
2556: }
2557: MatDenseRestoreArrayRead(A,&aa);
2558: VecRestoreArray(v,&x);
2559: return 0;
2560: }
2562: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2563: {
2564: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2565: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2566: PetscScalar *x;
2567: const PetscScalar *aa;
2570: MatDenseGetArrayRead(A,&aa);
2571: VecGetArray(v,&x);
2572: VecGetLocalSize(v,&p);
2574: for (i=0; i<m; i++) {
2575: x[i] = aa[i]; if (idx) idx[i] = 0;
2576: for (j=1; j<n; j++) {
2577: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2578: }
2579: }
2580: VecRestoreArray(v,&x);
2581: MatDenseRestoreArrayRead(A,&aa);
2582: return 0;
2583: }
2585: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2586: {
2587: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2588: PetscScalar *x;
2589: const PetscScalar *aa;
2592: MatDenseGetArrayRead(A,&aa);
2593: VecGetArray(v,&x);
2594: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2595: VecRestoreArray(v,&x);
2596: MatDenseRestoreArrayRead(A,&aa);
2597: return 0;
2598: }
2600: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions)
2601: {
2602: PetscInt i,j,m,n;
2603: const PetscScalar *a;
2605: MatGetSize(A,&m,&n);
2606: PetscArrayzero(reductions,n);
2607: MatDenseGetArrayRead(A,&a);
2608: if (type == NORM_2) {
2609: for (i=0; i<n; i++) {
2610: for (j=0; j<m; j++) {
2611: reductions[i] += PetscAbsScalar(a[j]*a[j]);
2612: }
2613: a += m;
2614: }
2615: } else if (type == NORM_1) {
2616: for (i=0; i<n; i++) {
2617: for (j=0; j<m; j++) {
2618: reductions[i] += PetscAbsScalar(a[j]);
2619: }
2620: a += m;
2621: }
2622: } else if (type == NORM_INFINITY) {
2623: for (i=0; i<n; i++) {
2624: for (j=0; j<m; j++) {
2625: reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]);
2626: }
2627: a += m;
2628: }
2629: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2630: for (i=0; i<n; i++) {
2631: for (j=0; j<m; j++) {
2632: reductions[i] += PetscRealPart(a[j]);
2633: }
2634: a += m;
2635: }
2636: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2637: for (i=0; i<n; i++) {
2638: for (j=0; j<m; j++) {
2639: reductions[i] += PetscImaginaryPart(a[j]);
2640: }
2641: a += m;
2642: }
2643: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type");
2644: MatDenseRestoreArrayRead(A,&a);
2645: if (type == NORM_2) {
2646: for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2647: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2648: for (i=0; i<n; i++) reductions[i] /= m;
2649: }
2650: return 0;
2651: }
2653: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2654: {
2655: PetscScalar *a;
2656: PetscInt lda,m,n,i,j;
2658: MatGetSize(x,&m,&n);
2659: MatDenseGetLDA(x,&lda);
2660: MatDenseGetArray(x,&a);
2661: for (j=0; j<n; j++) {
2662: for (i=0; i<m; i++) {
2663: PetscRandomGetValue(rctx,a+j*lda+i);
2664: }
2665: }
2666: MatDenseRestoreArray(x,&a);
2667: return 0;
2668: }
2670: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2671: {
2672: *missing = PETSC_FALSE;
2673: return 0;
2674: }
2676: /* vals is not const */
2677: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2678: {
2679: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2680: PetscScalar *v;
2683: MatDenseGetArray(A,&v);
2684: *vals = v+col*a->lda;
2685: MatDenseRestoreArray(A,&v);
2686: return 0;
2687: }
2689: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2690: {
2691: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
2692: return 0;
2693: }
2695: /* -------------------------------------------------------------------*/
2696: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2697: MatGetRow_SeqDense,
2698: MatRestoreRow_SeqDense,
2699: MatMult_SeqDense,
2700: /* 4*/ MatMultAdd_SeqDense,
2701: MatMultTranspose_SeqDense,
2702: MatMultTransposeAdd_SeqDense,
2703: NULL,
2704: NULL,
2705: NULL,
2706: /* 10*/ NULL,
2707: MatLUFactor_SeqDense,
2708: MatCholeskyFactor_SeqDense,
2709: MatSOR_SeqDense,
2710: MatTranspose_SeqDense,
2711: /* 15*/ MatGetInfo_SeqDense,
2712: MatEqual_SeqDense,
2713: MatGetDiagonal_SeqDense,
2714: MatDiagonalScale_SeqDense,
2715: MatNorm_SeqDense,
2716: /* 20*/ MatAssemblyBegin_SeqDense,
2717: MatAssemblyEnd_SeqDense,
2718: MatSetOption_SeqDense,
2719: MatZeroEntries_SeqDense,
2720: /* 24*/ MatZeroRows_SeqDense,
2721: NULL,
2722: NULL,
2723: NULL,
2724: NULL,
2725: /* 29*/ MatSetUp_SeqDense,
2726: NULL,
2727: NULL,
2728: NULL,
2729: NULL,
2730: /* 34*/ MatDuplicate_SeqDense,
2731: NULL,
2732: NULL,
2733: NULL,
2734: NULL,
2735: /* 39*/ MatAXPY_SeqDense,
2736: MatCreateSubMatrices_SeqDense,
2737: NULL,
2738: MatGetValues_SeqDense,
2739: MatCopy_SeqDense,
2740: /* 44*/ MatGetRowMax_SeqDense,
2741: MatScale_SeqDense,
2742: MatShift_Basic,
2743: NULL,
2744: MatZeroRowsColumns_SeqDense,
2745: /* 49*/ MatSetRandom_SeqDense,
2746: NULL,
2747: NULL,
2748: NULL,
2749: NULL,
2750: /* 54*/ NULL,
2751: NULL,
2752: NULL,
2753: NULL,
2754: NULL,
2755: /* 59*/ MatCreateSubMatrix_SeqDense,
2756: MatDestroy_SeqDense,
2757: MatView_SeqDense,
2758: NULL,
2759: NULL,
2760: /* 64*/ NULL,
2761: NULL,
2762: NULL,
2763: NULL,
2764: NULL,
2765: /* 69*/ MatGetRowMaxAbs_SeqDense,
2766: NULL,
2767: NULL,
2768: NULL,
2769: NULL,
2770: /* 74*/ NULL,
2771: NULL,
2772: NULL,
2773: NULL,
2774: NULL,
2775: /* 79*/ NULL,
2776: NULL,
2777: NULL,
2778: NULL,
2779: /* 83*/ MatLoad_SeqDense,
2780: MatIsSymmetric_SeqDense,
2781: MatIsHermitian_SeqDense,
2782: NULL,
2783: NULL,
2784: NULL,
2785: /* 89*/ NULL,
2786: NULL,
2787: MatMatMultNumeric_SeqDense_SeqDense,
2788: NULL,
2789: NULL,
2790: /* 94*/ NULL,
2791: NULL,
2792: NULL,
2793: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2794: NULL,
2795: /* 99*/ MatProductSetFromOptions_SeqDense,
2796: NULL,
2797: NULL,
2798: MatConjugate_SeqDense,
2799: NULL,
2800: /*104*/ NULL,
2801: MatRealPart_SeqDense,
2802: MatImaginaryPart_SeqDense,
2803: NULL,
2804: NULL,
2805: /*109*/ NULL,
2806: NULL,
2807: MatGetRowMin_SeqDense,
2808: MatGetColumnVector_SeqDense,
2809: MatMissingDiagonal_SeqDense,
2810: /*114*/ NULL,
2811: NULL,
2812: NULL,
2813: NULL,
2814: NULL,
2815: /*119*/ NULL,
2816: NULL,
2817: NULL,
2818: NULL,
2819: NULL,
2820: /*124*/ NULL,
2821: MatGetColumnReductions_SeqDense,
2822: NULL,
2823: NULL,
2824: NULL,
2825: /*129*/ NULL,
2826: NULL,
2827: NULL,
2828: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2829: NULL,
2830: /*134*/ NULL,
2831: NULL,
2832: NULL,
2833: NULL,
2834: NULL,
2835: /*139*/ NULL,
2836: NULL,
2837: NULL,
2838: NULL,
2839: NULL,
2840: MatCreateMPIMatConcatenateSeqMat_SeqDense,
2841: /*145*/ NULL,
2842: NULL,
2843: NULL
2844: };
2846: /*@C
2847: MatCreateSeqDense - Creates a sequential dense matrix that
2848: is stored in column major order (the usual Fortran 77 manner). Many
2849: of the matrix operations use the BLAS and LAPACK routines.
2851: Collective
2853: Input Parameters:
2854: + comm - MPI communicator, set to PETSC_COMM_SELF
2855: . m - number of rows
2856: . n - number of columns
2857: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2858: to control all matrix memory allocation.
2860: Output Parameter:
2861: . A - the matrix
2863: Notes:
2864: The data input variable is intended primarily for Fortran programmers
2865: who wish to allocate their own matrix memory space. Most users should
2866: set data=NULL.
2868: Level: intermediate
2870: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2871: @*/
2872: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2873: {
2874: MatCreate(comm,A);
2875: MatSetSizes(*A,m,n,m,n);
2876: MatSetType(*A,MATSEQDENSE);
2877: MatSeqDenseSetPreallocation(*A,data);
2878: return 0;
2879: }
2881: /*@C
2882: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2884: Collective
2886: Input Parameters:
2887: + B - the matrix
2888: - data - the array (or NULL)
2890: Notes:
2891: The data input variable is intended primarily for Fortran programmers
2892: who wish to allocate their own matrix memory space. Most users should
2893: need not call this routine.
2895: Level: intermediate
2897: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
2899: @*/
2900: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2901: {
2903: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2904: return 0;
2905: }
2907: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2908: {
2909: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2912: B->preallocated = PETSC_TRUE;
2914: PetscLayoutSetUp(B->rmap);
2915: PetscLayoutSetUp(B->cmap);
2917: if (b->lda <= 0) b->lda = B->rmap->n;
2919: if (!data) { /* petsc-allocated storage */
2920: if (!b->user_alloc) PetscFree(b->v);
2921: PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
2922: PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));
2924: b->user_alloc = PETSC_FALSE;
2925: } else { /* user-allocated storage */
2926: if (!b->user_alloc) PetscFree(b->v);
2927: b->v = data;
2928: b->user_alloc = PETSC_TRUE;
2929: }
2930: B->assembled = PETSC_TRUE;
2931: return 0;
2932: }
2934: #if defined(PETSC_HAVE_ELEMENTAL)
2935: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2936: {
2937: Mat mat_elemental;
2938: const PetscScalar *array;
2939: PetscScalar *v_colwise;
2940: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2942: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2943: MatDenseGetArrayRead(A,&array);
2944: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2945: k = 0;
2946: for (j=0; j<N; j++) {
2947: cols[j] = j;
2948: for (i=0; i<M; i++) {
2949: v_colwise[j*M+i] = array[k++];
2950: }
2951: }
2952: for (i=0; i<M; i++) {
2953: rows[i] = i;
2954: }
2955: MatDenseRestoreArrayRead(A,&array);
2957: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2958: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2959: MatSetType(mat_elemental,MATELEMENTAL);
2960: MatSetUp(mat_elemental);
2962: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2963: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2964: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2965: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2966: PetscFree3(v_colwise,rows,cols);
2968: if (reuse == MAT_INPLACE_MATRIX) {
2969: MatHeaderReplace(A,&mat_elemental);
2970: } else {
2971: *newmat = mat_elemental;
2972: }
2973: return 0;
2974: }
2975: #endif
2977: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2978: {
2979: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2980: PetscBool data;
2982: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2985: b->lda = lda;
2986: return 0;
2987: }
2989: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2990: {
2991: PetscMPIInt size;
2993: MPI_Comm_size(comm,&size);
2994: if (size == 1) {
2995: if (scall == MAT_INITIAL_MATRIX) {
2996: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2997: } else {
2998: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2999: }
3000: } else {
3001: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
3002: }
3003: return 0;
3004: }
3006: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3007: {
3008: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3012: if (!a->cvec) {
3013: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3014: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3015: }
3016: a->vecinuse = col + 1;
3017: MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3018: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3019: *v = a->cvec;
3020: return 0;
3021: }
3023: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3024: {
3025: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3029: a->vecinuse = 0;
3030: MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3031: VecResetArray(a->cvec);
3032: if (v) *v = NULL;
3033: return 0;
3034: }
3036: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3037: {
3038: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3042: if (!a->cvec) {
3043: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3044: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3045: }
3046: a->vecinuse = col + 1;
3047: MatDenseGetArrayRead(A,&a->ptrinuse);
3048: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3049: VecLockReadPush(a->cvec);
3050: *v = a->cvec;
3051: return 0;
3052: }
3054: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3055: {
3056: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3060: a->vecinuse = 0;
3061: MatDenseRestoreArrayRead(A,&a->ptrinuse);
3062: VecLockReadPop(a->cvec);
3063: VecResetArray(a->cvec);
3064: if (v) *v = NULL;
3065: return 0;
3066: }
3068: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3069: {
3070: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3074: if (!a->cvec) {
3075: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3076: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3077: }
3078: a->vecinuse = col + 1;
3079: MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3080: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3081: *v = a->cvec;
3082: return 0;
3083: }
3085: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3086: {
3087: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3091: a->vecinuse = 0;
3092: MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3093: VecResetArray(a->cvec);
3094: if (v) *v = NULL;
3095: return 0;
3096: }
3098: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3099: {
3100: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3104: if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3105: MatDestroy(&a->cmat);
3106: }
3107: if (!a->cmat) {
3108: MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,a->v+(size_t)cbegin*a->lda,&a->cmat);
3109: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3110: } else {
3111: MatDensePlaceArray(a->cmat,a->v+(size_t)cbegin*a->lda);
3112: }
3113: MatDenseSetLDA(a->cmat,a->lda);
3114: a->matinuse = cbegin + 1;
3115: *v = a->cmat;
3116: #if defined(PETSC_HAVE_CUDA)
3117: A->offloadmask = PETSC_OFFLOAD_CPU;
3118: #endif
3119: return 0;
3120: }
3122: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3123: {
3124: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3129: a->matinuse = 0;
3130: MatDenseResetArray(a->cmat);
3131: if (v) *v = NULL;
3132: return 0;
3133: }
3135: /*MC
3136: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3138: Options Database Keys:
3139: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3141: Level: beginner
3143: .seealso: MatCreateSeqDense()
3145: M*/
3146: PetscErrorCode MatCreate_SeqDense(Mat B)
3147: {
3148: Mat_SeqDense *b;
3149: PetscMPIInt size;
3151: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3154: PetscNewLog(B,&b);
3155: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3156: B->data = (void*)b;
3158: b->roworiented = PETSC_TRUE;
3160: PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3161: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3162: PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3163: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3164: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3165: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3166: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3167: PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3168: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3169: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3170: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3171: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3172: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3173: #if defined(PETSC_HAVE_ELEMENTAL)
3174: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3175: #endif
3176: #if defined(PETSC_HAVE_SCALAPACK)
3177: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3178: #endif
3179: #if defined(PETSC_HAVE_CUDA)
3180: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3181: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3182: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3183: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3184: #endif
3185: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3186: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3187: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3188: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3189: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3191: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3192: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3193: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3194: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3195: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3196: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3197: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3198: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3199: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3200: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3201: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3202: return 0;
3203: }
3205: /*@C
3206: 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.
3208: Not Collective
3210: Input Parameters:
3211: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3212: - col - column index
3214: Output Parameter:
3215: . vals - pointer to the data
3217: Level: intermediate
3219: .seealso: MatDenseRestoreColumn()
3220: @*/
3221: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3222: {
3226: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3227: return 0;
3228: }
3230: /*@C
3231: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3233: Not Collective
3235: Input Parameters:
3236: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3237: - vals - pointer to the data (may be NULL)
3239: Level: intermediate
3241: .seealso: MatDenseGetColumn()
3242: @*/
3243: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3244: {
3247: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3248: return 0;
3249: }
3251: /*@
3252: MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3254: Collective
3256: Input Parameters:
3257: + mat - the Mat object
3258: - col - the column index
3260: Output Parameter:
3261: . v - the vector
3263: Notes:
3264: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3265: Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.
3267: Level: intermediate
3269: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3270: @*/
3271: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3272: {
3279: PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3280: return 0;
3281: }
3283: /*@
3284: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3286: Collective
3288: Input Parameters:
3289: + mat - the Mat object
3290: . col - the column index
3291: - v - the Vec object (may be NULL)
3293: Level: intermediate
3295: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3296: @*/
3297: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3298: {
3304: PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3305: return 0;
3306: }
3308: /*@
3309: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3311: Collective
3313: Input Parameters:
3314: + mat - the Mat object
3315: - col - the column index
3317: Output Parameter:
3318: . v - the vector
3320: Notes:
3321: The vector is owned by PETSc and users cannot modify it.
3322: Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3323: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.
3325: Level: intermediate
3327: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3328: @*/
3329: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3330: {
3337: PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3338: return 0;
3339: }
3341: /*@
3342: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3344: Collective
3346: Input Parameters:
3347: + mat - the Mat object
3348: . col - the column index
3349: - v - the Vec object (may be NULL)
3351: Level: intermediate
3353: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3354: @*/
3355: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3356: {
3362: PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3363: return 0;
3364: }
3366: /*@
3367: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3369: Collective
3371: Input Parameters:
3372: + mat - the Mat object
3373: - col - the column index
3375: Output Parameter:
3376: . v - the vector
3378: Notes:
3379: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3380: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.
3382: Level: intermediate
3384: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3385: @*/
3386: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3387: {
3394: PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3395: return 0;
3396: }
3398: /*@
3399: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3401: Collective
3403: Input Parameters:
3404: + mat - the Mat object
3405: . col - the column index
3406: - v - the Vec object (may be NULL)
3408: Level: intermediate
3410: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3411: @*/
3412: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3413: {
3419: PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3420: return 0;
3421: }
3423: /*@
3424: MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3426: Collective
3428: Input Parameters:
3429: + mat - the Mat object
3430: . cbegin - the first index in the block
3431: - cend - the last index in the block
3433: Output Parameter:
3434: . v - the matrix
3436: Notes:
3437: The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3439: Level: intermediate
3441: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3442: @*/
3443: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3444: {
3453: PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3454: return 0;
3455: }
3457: /*@
3458: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3460: Collective
3462: Input Parameters:
3463: + mat - the Mat object
3464: - v - the Mat object (may be NULL)
3466: Level: intermediate
3468: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3469: @*/
3470: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3471: {
3475: PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3476: return 0;
3477: }