Actual source code: dense.c
petsc-3.12.5 2020-03-29
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;
19: if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
20: MatDenseGetArray(A,&v);
21: if (!hermitian) {
22: for (k=0;k<n;k++) {
23: for (j=k;j<n;j++) {
24: v[j*mat->lda + k] = v[k*mat->lda + j];
25: }
26: }
27: } else {
28: for (k=0;k<n;k++) {
29: for (j=k;j<n;j++) {
30: v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]);
31: }
32: }
33: }
34: MatDenseRestoreArray(A,&v);
35: return(0);
36: }
38: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
39: {
40: #if defined(PETSC_MISSING_LAPACK_POTRF)
42: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
43: #else
44: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
46: PetscBLASInt info,n;
49: if (!A->rmap->n || !A->cmap->n) return(0);
50: PetscBLASIntCast(A->cmap->n,&n);
51: if (A->factortype == MAT_FACTOR_LU) {
52: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
53: if (!mat->fwork) {
54: mat->lfwork = n;
55: PetscMalloc1(mat->lfwork,&mat->fwork);
56: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
57: }
58: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
59: PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
60: PetscFPTrapPop();
61: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
62: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
63: if (A->spd) {
64: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
65: PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
66: PetscFPTrapPop();
67: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
68: #if defined(PETSC_USE_COMPLEX)
69: } else if (A->hermitian) {
70: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
71: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
72: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
73: PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
74: PetscFPTrapPop();
75: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
76: #endif
77: } else { /* symmetric case */
78: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
79: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
80: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
81: #if defined(PETSC_MISSING_LAPACK_SYTRI)
82: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"SYTRI - Lapack routine is unavailable.");
83: #else
84: PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
85: #endif
86: PetscFPTrapPop();
87: MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
88: }
89: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
90: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
91: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
92: #endif
94: A->ops->solve = NULL;
95: A->ops->matsolve = NULL;
96: A->ops->solvetranspose = NULL;
97: A->ops->matsolvetranspose = NULL;
98: A->ops->solveadd = NULL;
99: A->ops->solvetransposeadd = NULL;
100: A->factortype = MAT_FACTOR_NONE;
101: PetscFree(A->solvertype);
102: return(0);
103: }
105: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
106: {
107: PetscErrorCode ierr;
108: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
109: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
110: PetscScalar *slot,*bb,*v;
111: const PetscScalar *xx;
114: #if defined(PETSC_USE_DEBUG)
115: for (i=0; i<N; i++) {
116: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
117: if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
118: if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n);
119: }
120: #endif
121: if (!N) return(0);
123: /* fix right hand side if needed */
124: if (x && b) {
125: Vec xt;
127: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
128: VecDuplicate(x,&xt);
129: VecCopy(x,xt);
130: VecScale(xt,-1.0);
131: MatMultAdd(A,xt,b,b);
132: VecDestroy(&xt);
133: VecGetArrayRead(x,&xx);
134: VecGetArray(b,&bb);
135: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
136: VecRestoreArrayRead(x,&xx);
137: VecRestoreArray(b,&bb);
138: }
140: MatDenseGetArray(A,&v);
141: for (i=0; i<N; i++) {
142: slot = v + rows[i]*m;
143: PetscArrayzero(slot,r);
144: }
145: for (i=0; i<N; i++) {
146: slot = v + rows[i];
147: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
148: }
149: if (diag != 0.0) {
150: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
151: for (i=0; i<N; i++) {
152: slot = v + (m+1)*rows[i];
153: *slot = diag;
154: }
155: }
156: MatDenseRestoreArray(A,&v);
157: return(0);
158: }
160: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
161: {
162: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
166: if (c->ptapwork) {
167: (*C->ops->matmultnumeric)(A,P,c->ptapwork);
168: (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
169: } else { /* first time went trough the basic. Should we add better dispatching for subclasses? */
170: MatPtAP_Basic(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
171: }
172: return(0);
173: }
175: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
176: {
177: Mat_SeqDense *c;
178: PetscBool flg;
182: PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
183: MatCreate(PetscObjectComm((PetscObject)A),C);
184: MatSetSizes(*C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
185: MatSetType(*C,flg ? ((PetscObject)A)->type_name : MATDENSE);
186: MatSeqDenseSetPreallocation(*C,NULL);
187: c = (Mat_SeqDense*)((*C)->data);
188: MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
189: MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
190: MatSetType(c->ptapwork,flg ? ((PetscObject)A)->type_name : MATDENSE);
191: MatSeqDenseSetPreallocation(c->ptapwork,NULL);
192: return(0);
193: }
195: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
196: {
200: if (reuse == MAT_INITIAL_MATRIX) {
201: MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
202: }
203: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
204: (*(*C)->ops->ptapnumeric)(A,P,*C);
205: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
206: return(0);
207: }
209: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
210: {
211: Mat B = NULL;
212: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
213: Mat_SeqDense *b;
215: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
216: MatScalar *av=a->a;
217: PetscBool isseqdense;
220: if (reuse == MAT_REUSE_MATRIX) {
221: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
222: if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
223: }
224: if (reuse != MAT_REUSE_MATRIX) {
225: MatCreate(PetscObjectComm((PetscObject)A),&B);
226: MatSetSizes(B,m,n,m,n);
227: MatSetType(B,MATSEQDENSE);
228: MatSeqDenseSetPreallocation(B,NULL);
229: b = (Mat_SeqDense*)(B->data);
230: } else {
231: b = (Mat_SeqDense*)((*newmat)->data);
232: PetscArrayzero(b->v,m*n);
233: }
234: for (i=0; i<m; i++) {
235: PetscInt j;
236: for (j=0;j<ai[1]-ai[0];j++) {
237: b->v[*aj*m+i] = *av;
238: aj++;
239: av++;
240: }
241: ai++;
242: }
244: if (reuse == MAT_INPLACE_MATRIX) {
245: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
246: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
247: MatHeaderReplace(A,&B);
248: } else {
249: if (B) *newmat = B;
250: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
251: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
252: }
253: return(0);
254: }
256: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
257: {
258: Mat B;
259: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
261: PetscInt i, j;
262: PetscInt *rows, *nnz;
263: MatScalar *aa = a->v, *vals;
266: MatCreate(PetscObjectComm((PetscObject)A),&B);
267: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
268: MatSetType(B,MATSEQAIJ);
269: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
270: for (j=0; j<A->cmap->n; j++) {
271: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
272: aa += a->lda;
273: }
274: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
275: aa = a->v;
276: for (j=0; j<A->cmap->n; j++) {
277: PetscInt numRows = 0;
278: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
279: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
280: aa += a->lda;
281: }
282: PetscFree3(rows,nnz,vals);
283: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
284: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
286: if (reuse == MAT_INPLACE_MATRIX) {
287: MatHeaderReplace(A,&B);
288: } else {
289: *newmat = B;
290: }
291: return(0);
292: }
294: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
295: {
296: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
297: const PetscScalar *xv;
298: PetscScalar *yv;
299: PetscBLASInt N,m,ldax,lday,one = 1;
300: PetscErrorCode ierr;
303: MatDenseGetArrayRead(X,&xv);
304: MatDenseGetArray(Y,&yv);
305: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
306: PetscBLASIntCast(X->rmap->n,&m);
307: PetscBLASIntCast(x->lda,&ldax);
308: PetscBLASIntCast(y->lda,&lday);
309: if (ldax>m || lday>m) {
310: PetscInt j;
312: for (j=0; j<X->cmap->n; j++) {
313: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
314: }
315: } else {
316: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
317: }
318: MatDenseRestoreArrayRead(X,&xv);
319: MatDenseRestoreArray(Y,&yv);
320: PetscLogFlops(PetscMax(2*N-1,0));
321: return(0);
322: }
324: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
325: {
326: PetscLogDouble N = A->rmap->n*A->cmap->n;
329: info->block_size = 1.0;
330: info->nz_allocated = N;
331: info->nz_used = N;
332: info->nz_unneeded = 0;
333: info->assemblies = A->num_ass;
334: info->mallocs = 0;
335: info->memory = ((PetscObject)A)->mem;
336: info->fill_ratio_given = 0;
337: info->fill_ratio_needed = 0;
338: info->factor_mallocs = 0;
339: return(0);
340: }
342: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
343: {
344: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
345: PetscScalar *v;
347: PetscBLASInt one = 1,j,nz,lda;
350: MatDenseGetArray(A,&v);
351: PetscBLASIntCast(a->lda,&lda);
352: if (lda>A->rmap->n) {
353: PetscBLASIntCast(A->rmap->n,&nz);
354: for (j=0; j<A->cmap->n; j++) {
355: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
356: }
357: } else {
358: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
359: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
360: }
361: PetscLogFlops(nz);
362: MatDenseRestoreArray(A,&v);
363: return(0);
364: }
366: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
367: {
368: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
369: PetscInt i,j,m = A->rmap->n,N = a->lda;
370: const PetscScalar *v;
371: PetscErrorCode ierr;
374: *fl = PETSC_FALSE;
375: if (A->rmap->n != A->cmap->n) return(0);
376: MatDenseGetArrayRead(A,&v);
377: for (i=0; i<m; i++) {
378: for (j=i; j<m; j++) {
379: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
380: }
381: }
382: MatDenseRestoreArrayRead(A,&v);
383: *fl = PETSC_TRUE;
384: return(0);
385: }
387: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
388: {
389: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
391: PetscInt lda = (PetscInt)mat->lda,j,m;
394: PetscLayoutReference(A->rmap,&newi->rmap);
395: PetscLayoutReference(A->cmap,&newi->cmap);
396: MatSeqDenseSetPreallocation(newi,NULL);
397: if (cpvalues == MAT_COPY_VALUES) {
398: const PetscScalar *av;
399: PetscScalar *v;
401: MatDenseGetArrayRead(A,&av);
402: MatDenseGetArray(newi,&v);
403: if (lda>A->rmap->n) {
404: m = A->rmap->n;
405: for (j=0; j<A->cmap->n; j++) {
406: PetscArraycpy(v+j*m,av+j*lda,m);
407: }
408: } else {
409: PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
410: }
411: MatDenseRestoreArray(newi,&v);
412: MatDenseRestoreArrayRead(A,&av);
413: }
414: return(0);
415: }
417: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
418: {
422: MatCreate(PetscObjectComm((PetscObject)A),newmat);
423: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
424: MatSetType(*newmat,((PetscObject)A)->type_name);
425: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
426: return(0);
427: }
429: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
430: {
431: MatFactorInfo info;
435: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
436: (*fact->ops->lufactor)(fact,0,0,&info);
437: return(0);
438: }
440: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
441: {
442: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
443: PetscErrorCode ierr;
444: const PetscScalar *x;
445: PetscScalar *y;
446: PetscBLASInt one = 1,info,m;
449: PetscBLASIntCast(A->rmap->n,&m);
450: VecGetArrayRead(xx,&x);
451: VecGetArray(yy,&y);
452: PetscArraycpy(y,x,A->rmap->n);
453: if (A->factortype == MAT_FACTOR_LU) {
454: #if defined(PETSC_MISSING_LAPACK_GETRS)
455: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
456: #else
457: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
458: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
459: PetscFPTrapPop();
460: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
461: #endif
462: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
463: #if defined(PETSC_MISSING_LAPACK_POTRS)
464: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
465: #else
466: if (A->spd) {
467: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
468: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
469: PetscFPTrapPop();
470: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
471: #if defined(PETSC_USE_COMPLEX)
472: } else if (A->hermitian) {
473: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
474: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
475: PetscFPTrapPop();
476: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
477: #endif
478: } else { /* symmetric case */
479: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
480: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
481: PetscFPTrapPop();
482: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
483: }
484: #endif
485: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
486: VecRestoreArrayRead(xx,&x);
487: VecRestoreArray(yy,&y);
488: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
489: return(0);
490: }
492: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
493: {
494: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
495: PetscErrorCode ierr;
496: const PetscScalar *b;
497: PetscScalar *x;
498: PetscInt n;
499: PetscBLASInt nrhs,info,m;
502: PetscBLASIntCast(A->rmap->n,&m);
503: MatGetSize(B,NULL,&n);
504: PetscBLASIntCast(n,&nrhs);
505: MatDenseGetArrayRead(B,&b);
506: MatDenseGetArray(X,&x);
508: PetscArraycpy(x,b,m*nrhs);
510: if (A->factortype == MAT_FACTOR_LU) {
511: #if defined(PETSC_MISSING_LAPACK_GETRS)
512: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
513: #else
514: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
515: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
516: PetscFPTrapPop();
517: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
518: #endif
519: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
520: #if defined(PETSC_MISSING_LAPACK_POTRS)
521: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
522: #else
523: if (A->spd) {
524: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
525: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
526: PetscFPTrapPop();
527: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
528: #if defined(PETSC_USE_COMPLEX)
529: } else if (A->hermitian) {
530: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
531: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
532: PetscFPTrapPop();
533: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
534: #endif
535: } else { /* symmetric case */
536: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
537: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
538: PetscFPTrapPop();
539: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
540: }
541: #endif
542: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
544: MatDenseRestoreArrayRead(B,&b);
545: MatDenseRestoreArray(X,&x);
546: PetscLogFlops(nrhs*(2.0*m*m - m));
547: return(0);
548: }
550: static PetscErrorCode MatConjugate_SeqDense(Mat);
552: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
553: {
554: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
555: PetscErrorCode ierr;
556: const PetscScalar *x;
557: PetscScalar *y;
558: PetscBLASInt one = 1,info,m;
561: PetscBLASIntCast(A->rmap->n,&m);
562: VecGetArrayRead(xx,&x);
563: VecGetArray(yy,&y);
564: PetscArraycpy(y,x,A->rmap->n);
565: if (A->factortype == MAT_FACTOR_LU) {
566: #if defined(PETSC_MISSING_LAPACK_GETRS)
567: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
568: #else
569: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
570: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
571: PetscFPTrapPop();
572: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
573: #endif
574: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
575: #if defined(PETSC_MISSING_LAPACK_POTRS)
576: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
577: #else
578: if (A->spd) {
579: #if defined(PETSC_USE_COMPLEX)
580: MatConjugate_SeqDense(A);
581: #endif
582: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
583: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
584: PetscFPTrapPop();
585: #if defined(PETSC_USE_COMPLEX)
586: MatConjugate_SeqDense(A);
587: #endif
588: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
589: #if defined(PETSC_USE_COMPLEX)
590: } else if (A->hermitian) {
591: MatConjugate_SeqDense(A);
592: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
593: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
594: PetscFPTrapPop();
595: MatConjugate_SeqDense(A);
596: #endif
597: } else { /* symmetric case */
598: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
599: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
600: PetscFPTrapPop();
601: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
602: }
603: #endif
604: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
605: VecRestoreArrayRead(xx,&x);
606: VecRestoreArray(yy,&y);
607: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
608: return(0);
609: }
611: /* ---------------------------------------------------------------*/
612: /* COMMENT: I have chosen to hide row permutation in the pivots,
613: rather than put it in the Mat->row slot.*/
614: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
615: {
616: #if defined(PETSC_MISSING_LAPACK_GETRF)
618: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
619: #else
620: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
622: PetscBLASInt n,m,info;
625: PetscBLASIntCast(A->cmap->n,&n);
626: PetscBLASIntCast(A->rmap->n,&m);
627: if (!mat->pivots) {
628: PetscMalloc1(A->rmap->n,&mat->pivots);
629: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
630: }
631: if (!A->rmap->n || !A->cmap->n) return(0);
632: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
633: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
634: PetscFPTrapPop();
636: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
637: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
639: A->ops->solve = MatSolve_SeqDense;
640: A->ops->matsolve = MatMatSolve_SeqDense;
641: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
642: A->factortype = MAT_FACTOR_LU;
644: PetscFree(A->solvertype);
645: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
647: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
648: #endif
649: return(0);
650: }
652: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
653: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
654: {
655: #if defined(PETSC_MISSING_LAPACK_POTRF)
657: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
658: #else
659: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
661: PetscBLASInt info,n;
664: PetscBLASIntCast(A->cmap->n,&n);
665: if (!A->rmap->n || !A->cmap->n) return(0);
666: if (A->spd) {
667: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
668: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
669: PetscFPTrapPop();
670: #if defined(PETSC_USE_COMPLEX)
671: } else if (A->hermitian) {
672: if (!mat->pivots) {
673: PetscMalloc1(A->rmap->n,&mat->pivots);
674: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
675: }
676: if (!mat->fwork) {
677: PetscScalar dummy;
679: mat->lfwork = -1;
680: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
681: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
682: PetscFPTrapPop();
683: mat->lfwork = (PetscInt)PetscRealPart(dummy);
684: PetscMalloc1(mat->lfwork,&mat->fwork);
685: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
686: }
687: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
688: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
689: PetscFPTrapPop();
690: #endif
691: } else { /* symmetric case */
692: if (!mat->pivots) {
693: PetscMalloc1(A->rmap->n,&mat->pivots);
694: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
695: }
696: if (!mat->fwork) {
697: PetscScalar dummy;
699: mat->lfwork = -1;
700: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
701: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
702: PetscFPTrapPop();
703: mat->lfwork = (PetscInt)PetscRealPart(dummy);
704: PetscMalloc1(mat->lfwork,&mat->fwork);
705: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
706: }
707: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
708: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
709: PetscFPTrapPop();
710: }
711: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
713: A->ops->solve = MatSolve_SeqDense;
714: A->ops->matsolve = MatMatSolve_SeqDense;
715: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
716: A->factortype = MAT_FACTOR_CHOLESKY;
718: PetscFree(A->solvertype);
719: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
721: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
722: #endif
723: return(0);
724: }
727: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
728: {
730: MatFactorInfo info;
733: info.fill = 1.0;
735: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
736: (*fact->ops->choleskyfactor)(fact,0,&info);
737: return(0);
738: }
740: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
741: {
743: fact->assembled = PETSC_TRUE;
744: fact->preallocated = PETSC_TRUE;
745: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
746: fact->ops->solve = MatSolve_SeqDense;
747: fact->ops->matsolve = MatMatSolve_SeqDense;
748: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
749: return(0);
750: }
752: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
753: {
755: fact->preallocated = PETSC_TRUE;
756: fact->assembled = PETSC_TRUE;
757: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
758: fact->ops->solve = MatSolve_SeqDense;
759: fact->ops->matsolve = MatMatSolve_SeqDense;
760: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
761: return(0);
762: }
764: /* uses LAPACK */
765: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
766: {
770: MatCreate(PetscObjectComm((PetscObject)A),fact);
771: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
772: MatSetType(*fact,MATDENSE);
773: if (ftype == MAT_FACTOR_LU) {
774: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
775: } else {
776: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
777: }
778: (*fact)->factortype = ftype;
780: PetscFree((*fact)->solvertype);
781: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
782: return(0);
783: }
785: /* ------------------------------------------------------------------*/
786: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
787: {
788: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
789: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
790: const PetscScalar *b;
791: PetscErrorCode ierr;
792: PetscInt m = A->rmap->n,i;
793: PetscBLASInt o = 1,bm;
796: #if defined(PETSC_HAVE_CUDA)
797: if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
798: #endif
799: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
800: PetscBLASIntCast(m,&bm);
801: if (flag & SOR_ZERO_INITIAL_GUESS) {
802: /* this is a hack fix, should have another version without the second BLASdotu */
803: VecSet(xx,zero);
804: }
805: VecGetArray(xx,&x);
806: VecGetArrayRead(bb,&b);
807: its = its*lits;
808: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
809: while (its--) {
810: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
811: for (i=0; i<m; i++) {
812: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
813: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
814: }
815: }
816: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
817: for (i=m-1; i>=0; i--) {
818: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
819: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
820: }
821: }
822: }
823: VecRestoreArrayRead(bb,&b);
824: VecRestoreArray(xx,&x);
825: return(0);
826: }
828: /* -----------------------------------------------------------------*/
829: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
830: {
831: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
832: const PetscScalar *v = mat->v,*x;
833: PetscScalar *y;
834: PetscErrorCode ierr;
835: PetscBLASInt m, n,_One=1;
836: PetscScalar _DOne=1.0,_DZero=0.0;
839: PetscBLASIntCast(A->rmap->n,&m);
840: PetscBLASIntCast(A->cmap->n,&n);
841: VecGetArrayRead(xx,&x);
842: VecGetArrayWrite(yy,&y);
843: if (!A->rmap->n || !A->cmap->n) {
844: PetscBLASInt i;
845: for (i=0; i<n; i++) y[i] = 0.0;
846: } else {
847: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
848: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
849: }
850: VecRestoreArrayRead(xx,&x);
851: VecRestoreArrayWrite(yy,&y);
852: return(0);
853: }
855: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
856: {
857: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
858: PetscScalar *y,_DOne=1.0,_DZero=0.0;
859: PetscErrorCode ierr;
860: PetscBLASInt m, n, _One=1;
861: const PetscScalar *v = mat->v,*x;
864: PetscBLASIntCast(A->rmap->n,&m);
865: PetscBLASIntCast(A->cmap->n,&n);
866: VecGetArrayRead(xx,&x);
867: VecGetArrayWrite(yy,&y);
868: if (!A->rmap->n || !A->cmap->n) {
869: PetscBLASInt i;
870: for (i=0; i<m; i++) y[i] = 0.0;
871: } else {
872: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
873: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
874: }
875: VecRestoreArrayRead(xx,&x);
876: VecRestoreArrayWrite(yy,&y);
877: return(0);
878: }
880: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
881: {
882: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
883: const PetscScalar *v = mat->v,*x;
884: PetscScalar *y,_DOne=1.0;
885: PetscErrorCode ierr;
886: PetscBLASInt m, n, _One=1;
889: PetscBLASIntCast(A->rmap->n,&m);
890: PetscBLASIntCast(A->cmap->n,&n);
891: if (!A->rmap->n || !A->cmap->n) return(0);
892: if (zz != yy) {VecCopy(zz,yy);}
893: VecGetArrayRead(xx,&x);
894: VecGetArray(yy,&y);
895: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
896: VecRestoreArrayRead(xx,&x);
897: VecRestoreArray(yy,&y);
898: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
899: return(0);
900: }
902: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
903: {
904: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
905: const PetscScalar *v = mat->v,*x;
906: PetscScalar *y;
907: PetscErrorCode ierr;
908: PetscBLASInt m, n, _One=1;
909: PetscScalar _DOne=1.0;
912: PetscBLASIntCast(A->rmap->n,&m);
913: PetscBLASIntCast(A->cmap->n,&n);
914: if (!A->rmap->n || !A->cmap->n) return(0);
915: if (zz != yy) {VecCopy(zz,yy);}
916: VecGetArrayRead(xx,&x);
917: VecGetArray(yy,&y);
918: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
919: VecRestoreArrayRead(xx,&x);
920: VecRestoreArray(yy,&y);
921: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
922: return(0);
923: }
925: /* -----------------------------------------------------------------*/
926: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
927: {
928: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
930: PetscInt i;
933: *ncols = A->cmap->n;
934: if (cols) {
935: PetscMalloc1(A->cmap->n+1,cols);
936: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
937: }
938: if (vals) {
939: const PetscScalar *v;
941: MatDenseGetArrayRead(A,&v);
942: PetscMalloc1(A->cmap->n+1,vals);
943: v += row;
944: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
945: MatDenseRestoreArrayRead(A,&v);
946: }
947: return(0);
948: }
950: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
951: {
955: if (cols) {PetscFree(*cols);}
956: if (vals) {PetscFree(*vals); }
957: return(0);
958: }
959: /* ----------------------------------------------------------------*/
960: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
961: {
962: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
963: PetscScalar *av;
964: PetscInt i,j,idx=0;
965: #if defined(PETSC_HAVE_CUDA)
966: PetscOffloadMask oldf;
967: #endif
968: PetscErrorCode ierr;
971: MatDenseGetArray(A,&av);
972: if (!mat->roworiented) {
973: if (addv == INSERT_VALUES) {
974: for (j=0; j<n; j++) {
975: if (indexn[j] < 0) {idx += m; continue;}
976: #if defined(PETSC_USE_DEBUG)
977: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
978: #endif
979: for (i=0; i<m; i++) {
980: if (indexm[i] < 0) {idx++; continue;}
981: #if defined(PETSC_USE_DEBUG)
982: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
983: #endif
984: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
985: }
986: }
987: } else {
988: for (j=0; j<n; j++) {
989: if (indexn[j] < 0) {idx += m; continue;}
990: #if defined(PETSC_USE_DEBUG)
991: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
992: #endif
993: for (i=0; i<m; i++) {
994: if (indexm[i] < 0) {idx++; continue;}
995: #if defined(PETSC_USE_DEBUG)
996: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
997: #endif
998: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
999: }
1000: }
1001: }
1002: } else {
1003: if (addv == INSERT_VALUES) {
1004: for (i=0; i<m; i++) {
1005: if (indexm[i] < 0) { idx += n; continue;}
1006: #if defined(PETSC_USE_DEBUG)
1007: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
1008: #endif
1009: for (j=0; j<n; j++) {
1010: if (indexn[j] < 0) { idx++; continue;}
1011: #if defined(PETSC_USE_DEBUG)
1012: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
1013: #endif
1014: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1015: }
1016: }
1017: } else {
1018: for (i=0; i<m; i++) {
1019: if (indexm[i] < 0) { idx += n; continue;}
1020: #if defined(PETSC_USE_DEBUG)
1021: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
1022: #endif
1023: for (j=0; j<n; j++) {
1024: if (indexn[j] < 0) { idx++; continue;}
1025: #if defined(PETSC_USE_DEBUG)
1026: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
1027: #endif
1028: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1029: }
1030: }
1031: }
1032: }
1033: /* hack to prevent unneeded copy to the GPU while returning the array */
1034: #if defined(PETSC_HAVE_CUDA)
1035: oldf = A->offloadmask;
1036: A->offloadmask = PETSC_OFFLOAD_GPU;
1037: #endif
1038: MatDenseRestoreArray(A,&av);
1039: #if defined(PETSC_HAVE_CUDA)
1040: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1041: #endif
1042: return(0);
1043: }
1045: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1046: {
1047: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1048: const PetscScalar *vv;
1049: PetscInt i,j;
1050: PetscErrorCode ierr;
1053: MatDenseGetArrayRead(A,&vv);
1054: /* row-oriented output */
1055: for (i=0; i<m; i++) {
1056: if (indexm[i] < 0) {v += n;continue;}
1057: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n);
1058: for (j=0; j<n; j++) {
1059: if (indexn[j] < 0) {v++; continue;}
1060: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n);
1061: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1062: }
1063: }
1064: MatDenseRestoreArrayRead(A,&vv);
1065: return(0);
1066: }
1068: /* -----------------------------------------------------------------*/
1070: static PetscErrorCode MatLoad_SeqDense_Binary(Mat newmat,PetscViewer viewer)
1071: {
1072: Mat_SeqDense *a;
1074: PetscInt *scols,i,j,nz,header[4];
1075: int fd;
1076: PetscMPIInt size;
1077: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
1078: PetscScalar *vals,*svals,*v,*w;
1079: MPI_Comm comm;
1082: PetscObjectGetComm((PetscObject)viewer,&comm);
1083: MPI_Comm_size(comm,&size);
1084: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1085: PetscViewerBinaryGetDescriptor(viewer,&fd);
1086: PetscBinaryRead(fd,header,4,NULL,PETSC_INT);
1087: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1088: M = header[1]; N = header[2]; nz = header[3];
1090: /* set global size if not set already*/
1091: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1092: MatSetSizes(newmat,M,N,M,N);
1093: } else {
1094: /* if sizes and type are already set, check if the vector global sizes are correct */
1095: MatGetSize(newmat,&grows,&gcols);
1096: if (M != grows || N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols);
1097: }
1098: a = (Mat_SeqDense*)newmat->data;
1099: if (!a->user_alloc) {
1100: MatSeqDenseSetPreallocation(newmat,NULL);
1101: }
1103: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1104: a = (Mat_SeqDense*)newmat->data;
1105: v = a->v;
1106: /* Allocate some temp space to read in the values and then flip them
1107: from row major to column major */
1108: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1109: /* read in nonzero values */
1110: PetscBinaryRead(fd,w,M*N,NULL,PETSC_SCALAR);
1111: /* now flip the values and store them in the matrix*/
1112: for (j=0; j<N; j++) {
1113: for (i=0; i<M; i++) {
1114: *v++ =w[i*N+j];
1115: }
1116: }
1117: PetscFree(w);
1118: } else {
1119: /* read row lengths */
1120: PetscMalloc1(M+1,&rowlengths);
1121: PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);
1123: a = (Mat_SeqDense*)newmat->data;
1124: v = a->v;
1126: /* read column indices and nonzeros */
1127: PetscMalloc1(nz+1,&scols);
1128: cols = scols;
1129: PetscBinaryRead(fd,cols,nz,NULL,PETSC_INT);
1130: PetscMalloc1(nz+1,&svals);
1131: vals = svals;
1132: PetscBinaryRead(fd,vals,nz,NULL,PETSC_SCALAR);
1134: /* insert into matrix */
1135: for (i=0; i<M; i++) {
1136: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1137: svals += rowlengths[i]; scols += rowlengths[i];
1138: }
1139: PetscFree(vals);
1140: PetscFree(cols);
1141: PetscFree(rowlengths);
1142: }
1143: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1144: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1145: return(0);
1146: }
1148: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1149: {
1150: PetscBool isbinary, ishdf5;
1156: /* force binary viewer to load .info file if it has not yet done so */
1157: PetscViewerSetUp(viewer);
1158: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1159: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1160: if (isbinary) {
1161: MatLoad_SeqDense_Binary(newMat,viewer);
1162: } else if (ishdf5) {
1163: #if defined(PETSC_HAVE_HDF5)
1164: MatLoad_Dense_HDF5(newMat,viewer);
1165: #else
1166: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1167: #endif
1168: } else {
1169: SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
1170: }
1171: return(0);
1172: }
1174: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1175: {
1176: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1177: PetscErrorCode ierr;
1178: PetscInt i,j;
1179: const char *name;
1180: PetscScalar *v,*av;
1181: PetscViewerFormat format;
1182: #if defined(PETSC_USE_COMPLEX)
1183: PetscBool allreal = PETSC_TRUE;
1184: #endif
1187: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1188: PetscViewerGetFormat(viewer,&format);
1189: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1190: return(0); /* do nothing for now */
1191: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1192: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1193: for (i=0; i<A->rmap->n; i++) {
1194: v = av + i;
1195: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1196: for (j=0; j<A->cmap->n; j++) {
1197: #if defined(PETSC_USE_COMPLEX)
1198: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1199: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1200: } else if (PetscRealPart(*v)) {
1201: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1202: }
1203: #else
1204: if (*v) {
1205: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1206: }
1207: #endif
1208: v += a->lda;
1209: }
1210: PetscViewerASCIIPrintf(viewer,"\n");
1211: }
1212: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1213: } else {
1214: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1215: #if defined(PETSC_USE_COMPLEX)
1216: /* determine if matrix has all real values */
1217: v = av;
1218: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1219: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1220: }
1221: #endif
1222: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1223: PetscObjectGetName((PetscObject)A,&name);
1224: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1225: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1226: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1227: }
1229: for (i=0; i<A->rmap->n; i++) {
1230: v = av + i;
1231: for (j=0; j<A->cmap->n; j++) {
1232: #if defined(PETSC_USE_COMPLEX)
1233: if (allreal) {
1234: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1235: } else {
1236: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1237: }
1238: #else
1239: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1240: #endif
1241: v += a->lda;
1242: }
1243: PetscViewerASCIIPrintf(viewer,"\n");
1244: }
1245: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1246: PetscViewerASCIIPrintf(viewer,"];\n");
1247: }
1248: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1249: }
1250: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1251: PetscViewerFlush(viewer);
1252: return(0);
1253: }
1255: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1256: {
1257: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1258: PetscErrorCode ierr;
1259: int fd;
1260: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1261: PetscScalar *av,*v,*anonz,*vals;
1262: PetscViewerFormat format;
1265: PetscViewerBinaryGetDescriptor(viewer,&fd);
1266: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1267: PetscViewerGetFormat(viewer,&format);
1268: if (format == PETSC_VIEWER_NATIVE) {
1269: /* store the matrix as a dense matrix */
1270: PetscMalloc1(4,&col_lens);
1272: col_lens[0] = MAT_FILE_CLASSID;
1273: col_lens[1] = m;
1274: col_lens[2] = n;
1275: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1277: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1278: PetscFree(col_lens);
1280: /* write out matrix, by rows */
1281: PetscMalloc1(m*n+1,&vals);
1282: v = av;
1283: for (j=0; j<n; j++) {
1284: for (i=0; i<m; i++) {
1285: vals[j + i*n] = *v++;
1286: }
1287: }
1288: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1289: PetscFree(vals);
1290: } else {
1291: PetscMalloc1(4+nz,&col_lens);
1293: col_lens[0] = MAT_FILE_CLASSID;
1294: col_lens[1] = m;
1295: col_lens[2] = n;
1296: col_lens[3] = nz;
1298: /* store lengths of each row and write (including header) to file */
1299: for (i=0; i<m; i++) col_lens[4+i] = n;
1300: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1302: /* Possibly should write in smaller increments, not whole matrix at once? */
1303: /* store column indices (zero start index) */
1304: ict = 0;
1305: for (i=0; i<m; i++) {
1306: for (j=0; j<n; j++) col_lens[ict++] = j;
1307: }
1308: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1309: PetscFree(col_lens);
1311: /* store nonzero values */
1312: PetscMalloc1(nz+1,&anonz);
1313: ict = 0;
1314: for (i=0; i<m; i++) {
1315: v = av + i;
1316: for (j=0; j<n; j++) {
1317: anonz[ict++] = *v; v += a->lda;
1318: }
1319: }
1320: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1321: PetscFree(anonz);
1322: }
1323: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1324: return(0);
1325: }
1327: #include <petscdraw.h>
1328: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1329: {
1330: Mat A = (Mat) Aa;
1331: PetscErrorCode ierr;
1332: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1333: int color = PETSC_DRAW_WHITE;
1334: const PetscScalar *v;
1335: PetscViewer viewer;
1336: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1337: PetscViewerFormat format;
1340: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1341: PetscViewerGetFormat(viewer,&format);
1342: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1344: /* Loop over matrix elements drawing boxes */
1345: MatDenseGetArrayRead(A,&v);
1346: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1347: PetscDrawCollectiveBegin(draw);
1348: /* Blue for negative and Red for positive */
1349: for (j = 0; j < n; j++) {
1350: x_l = j; x_r = x_l + 1.0;
1351: for (i = 0; i < m; i++) {
1352: y_l = m - i - 1.0;
1353: y_r = y_l + 1.0;
1354: if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED;
1355: else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE;
1356: else continue;
1357: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1358: }
1359: }
1360: PetscDrawCollectiveEnd(draw);
1361: } else {
1362: /* use contour shading to indicate magnitude of values */
1363: /* first determine max of all nonzero values */
1364: PetscReal minv = 0.0, maxv = 0.0;
1365: PetscDraw popup;
1367: for (i=0; i < m*n; i++) {
1368: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1369: }
1370: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1371: PetscDrawGetPopup(draw,&popup);
1372: PetscDrawScalePopup(popup,minv,maxv);
1374: PetscDrawCollectiveBegin(draw);
1375: for (j=0; j<n; j++) {
1376: x_l = j;
1377: x_r = x_l + 1.0;
1378: for (i=0; i<m; i++) {
1379: y_l = m - i - 1.0;
1380: y_r = y_l + 1.0;
1381: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1382: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1383: }
1384: }
1385: PetscDrawCollectiveEnd(draw);
1386: }
1387: MatDenseRestoreArrayRead(A,&v);
1388: return(0);
1389: }
1391: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1392: {
1393: PetscDraw draw;
1394: PetscBool isnull;
1395: PetscReal xr,yr,xl,yl,h,w;
1399: PetscViewerDrawGetDraw(viewer,0,&draw);
1400: PetscDrawIsNull(draw,&isnull);
1401: if (isnull) return(0);
1403: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1404: xr += w; yr += h; xl = -w; yl = -h;
1405: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1406: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1407: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1408: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1409: PetscDrawSave(draw);
1410: return(0);
1411: }
1413: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1414: {
1416: PetscBool iascii,isbinary,isdraw;
1419: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1420: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1421: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1423: if (iascii) {
1424: MatView_SeqDense_ASCII(A,viewer);
1425: } else if (isbinary) {
1426: MatView_SeqDense_Binary(A,viewer);
1427: } else if (isdraw) {
1428: MatView_SeqDense_Draw(A,viewer);
1429: }
1430: return(0);
1431: }
1433: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1434: {
1435: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1438: a->unplacedarray = a->v;
1439: a->unplaced_user_alloc = a->user_alloc;
1440: a->v = (PetscScalar*) array;
1441: #if defined(PETSC_HAVE_CUDA)
1442: A->offloadmask = PETSC_OFFLOAD_CPU;
1443: #endif
1444: return(0);
1445: }
1447: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1448: {
1449: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1452: a->v = a->unplacedarray;
1453: a->user_alloc = a->unplaced_user_alloc;
1454: a->unplacedarray = NULL;
1455: #if defined(PETSC_HAVE_CUDA)
1456: A->offloadmask = PETSC_OFFLOAD_CPU;
1457: #endif
1458: return(0);
1459: }
1461: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1462: {
1463: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1467: #if defined(PETSC_USE_LOG)
1468: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1469: #endif
1470: PetscFree(l->pivots);
1471: PetscFree(l->fwork);
1472: MatDestroy(&l->ptapwork);
1473: if (!l->user_alloc) {PetscFree(l->v);}
1474: PetscFree(mat->data);
1476: PetscObjectChangeTypeName((PetscObject)mat,0);
1477: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1478: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1479: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1480: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1481: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1482: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1483: #if defined(PETSC_HAVE_ELEMENTAL)
1484: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1485: #endif
1486: #if defined(PETSC_HAVE_CUDA)
1487: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1488: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijcusparse_seqdense_C",NULL);
1489: #endif
1490: #if defined(PETSC_HAVE_VIENNACL)
1491: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijviennacl_seqdense_C",NULL);
1492: #endif
1493: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1494: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1495: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1496: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1497: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqbaij_seqdense_C",NULL);
1498: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);
1499: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);
1500: PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1501: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1502: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1503: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1504: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1505: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1506: return(0);
1507: }
1509: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1510: {
1511: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1513: PetscInt k,j,m,n,M;
1514: PetscScalar *v,tmp;
1517: m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1518: if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1519: MatDenseGetArray(A,&v);
1520: for (j=0; j<m; j++) {
1521: for (k=0; k<j; k++) {
1522: tmp = v[j + k*M];
1523: v[j + k*M] = v[k + j*M];
1524: v[k + j*M] = tmp;
1525: }
1526: }
1527: MatDenseRestoreArray(A,&v);
1528: } else { /* out-of-place transpose */
1529: Mat tmat;
1530: Mat_SeqDense *tmatd;
1531: PetscScalar *v2;
1532: PetscInt M2;
1534: if (reuse != MAT_REUSE_MATRIX) {
1535: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1536: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1537: MatSetType(tmat,((PetscObject)A)->type_name);
1538: MatSeqDenseSetPreallocation(tmat,NULL);
1539: } else tmat = *matout;
1541: MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1542: MatDenseGetArray(tmat,&v2);
1543: tmatd = (Mat_SeqDense*)tmat->data;
1544: M2 = tmatd->lda;
1545: for (j=0; j<n; j++) {
1546: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1547: }
1548: MatDenseRestoreArray(tmat,&v2);
1549: MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1550: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1551: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1552: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1553: else {
1554: MatHeaderMerge(A,&tmat);
1555: }
1556: }
1557: return(0);
1558: }
1560: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1561: {
1562: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1563: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1564: PetscInt i;
1565: const PetscScalar *v1,*v2;
1566: PetscErrorCode ierr;
1569: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1570: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1571: MatDenseGetArrayRead(A1,&v1);
1572: MatDenseGetArrayRead(A2,&v2);
1573: for (i=0; i<A1->cmap->n; i++) {
1574: PetscArraycmp(v1,v2,A1->rmap->n,flg);
1575: if (*flg == PETSC_FALSE) return(0);
1576: v1 += mat1->lda;
1577: v2 += mat2->lda;
1578: }
1579: MatDenseRestoreArrayRead(A1,&v1);
1580: MatDenseRestoreArrayRead(A2,&v2);
1581: *flg = PETSC_TRUE;
1582: return(0);
1583: }
1585: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1586: {
1587: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1588: PetscInt i,n,len;
1589: PetscScalar *x;
1590: const PetscScalar *vv;
1591: PetscErrorCode ierr;
1594: VecGetSize(v,&n);
1595: VecGetArray(v,&x);
1596: len = PetscMin(A->rmap->n,A->cmap->n);
1597: MatDenseGetArrayRead(A,&vv);
1598: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1599: for (i=0; i<len; i++) {
1600: x[i] = vv[i*mat->lda + i];
1601: }
1602: MatDenseRestoreArrayRead(A,&vv);
1603: VecRestoreArray(v,&x);
1604: return(0);
1605: }
1607: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1608: {
1609: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1610: const PetscScalar *l,*r;
1611: PetscScalar x,*v,*vv;
1612: PetscErrorCode ierr;
1613: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1616: MatDenseGetArray(A,&vv);
1617: if (ll) {
1618: VecGetSize(ll,&m);
1619: VecGetArrayRead(ll,&l);
1620: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1621: for (i=0; i<m; i++) {
1622: x = l[i];
1623: v = vv + i;
1624: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1625: }
1626: VecRestoreArrayRead(ll,&l);
1627: PetscLogFlops(1.0*n*m);
1628: }
1629: if (rr) {
1630: VecGetSize(rr,&n);
1631: VecGetArrayRead(rr,&r);
1632: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1633: for (i=0; i<n; i++) {
1634: x = r[i];
1635: v = vv + i*mat->lda;
1636: for (j=0; j<m; j++) (*v++) *= x;
1637: }
1638: VecRestoreArrayRead(rr,&r);
1639: PetscLogFlops(1.0*n*m);
1640: }
1641: MatDenseRestoreArray(A,&vv);
1642: return(0);
1643: }
1645: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1646: {
1647: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1648: PetscScalar *v,*vv;
1649: PetscReal sum = 0.0;
1650: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1651: PetscErrorCode ierr;
1654: MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1655: v = vv;
1656: if (type == NORM_FROBENIUS) {
1657: if (lda>m) {
1658: for (j=0; j<A->cmap->n; j++) {
1659: v = vv+j*lda;
1660: for (i=0; i<m; i++) {
1661: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1662: }
1663: }
1664: } else {
1665: #if defined(PETSC_USE_REAL___FP16)
1666: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1667: *nrm = BLASnrm2_(&cnt,v,&one);
1668: }
1669: #else
1670: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1671: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1672: }
1673: }
1674: *nrm = PetscSqrtReal(sum);
1675: #endif
1676: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1677: } else if (type == NORM_1) {
1678: *nrm = 0.0;
1679: for (j=0; j<A->cmap->n; j++) {
1680: v = vv + j*mat->lda;
1681: sum = 0.0;
1682: for (i=0; i<A->rmap->n; i++) {
1683: sum += PetscAbsScalar(*v); v++;
1684: }
1685: if (sum > *nrm) *nrm = sum;
1686: }
1687: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1688: } else if (type == NORM_INFINITY) {
1689: *nrm = 0.0;
1690: for (j=0; j<A->rmap->n; j++) {
1691: v = vv + j;
1692: sum = 0.0;
1693: for (i=0; i<A->cmap->n; i++) {
1694: sum += PetscAbsScalar(*v); v += mat->lda;
1695: }
1696: if (sum > *nrm) *nrm = sum;
1697: }
1698: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1699: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1700: MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1701: return(0);
1702: }
1704: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1705: {
1706: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1710: switch (op) {
1711: case MAT_ROW_ORIENTED:
1712: aij->roworiented = flg;
1713: break;
1714: case MAT_NEW_NONZERO_LOCATIONS:
1715: case MAT_NEW_NONZERO_LOCATION_ERR:
1716: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1717: case MAT_NEW_DIAGONALS:
1718: case MAT_KEEP_NONZERO_PATTERN:
1719: case MAT_IGNORE_OFF_PROC_ENTRIES:
1720: case MAT_USE_HASH_TABLE:
1721: case MAT_IGNORE_ZERO_ENTRIES:
1722: case MAT_IGNORE_LOWER_TRIANGULAR:
1723: case MAT_SORTED_FULL:
1724: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1725: break;
1726: case MAT_SPD:
1727: case MAT_SYMMETRIC:
1728: case MAT_STRUCTURALLY_SYMMETRIC:
1729: case MAT_HERMITIAN:
1730: case MAT_SYMMETRY_ETERNAL:
1731: /* These options are handled directly by MatSetOption() */
1732: break;
1733: default:
1734: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1735: }
1736: return(0);
1737: }
1739: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1740: {
1741: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1743: PetscInt lda=l->lda,m=A->rmap->n,j;
1744: PetscScalar *v;
1747: MatDenseGetArray(A,&v);
1748: if (lda>m) {
1749: for (j=0; j<A->cmap->n; j++) {
1750: PetscArrayzero(v+j*lda,m);
1751: }
1752: } else {
1753: PetscArrayzero(v,A->rmap->n*A->cmap->n);
1754: }
1755: MatDenseRestoreArray(A,&v);
1756: return(0);
1757: }
1759: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1760: {
1761: PetscErrorCode ierr;
1762: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1763: PetscInt m = l->lda, n = A->cmap->n, i,j;
1764: PetscScalar *slot,*bb,*v;
1765: const PetscScalar *xx;
1768: #if defined(PETSC_USE_DEBUG)
1769: for (i=0; i<N; i++) {
1770: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1771: if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
1772: }
1773: #endif
1774: if (!N) return(0);
1776: /* fix right hand side if needed */
1777: if (x && b) {
1778: VecGetArrayRead(x,&xx);
1779: VecGetArray(b,&bb);
1780: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1781: VecRestoreArrayRead(x,&xx);
1782: VecRestoreArray(b,&bb);
1783: }
1785: MatDenseGetArray(A,&v);
1786: for (i=0; i<N; i++) {
1787: slot = v + rows[i];
1788: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1789: }
1790: if (diag != 0.0) {
1791: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1792: for (i=0; i<N; i++) {
1793: slot = v + (m+1)*rows[i];
1794: *slot = diag;
1795: }
1796: }
1797: MatDenseRestoreArray(A,&v);
1798: return(0);
1799: }
1801: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1802: {
1803: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1806: *lda = mat->lda;
1807: return(0);
1808: }
1810: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1811: {
1812: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1815: *array = mat->v;
1816: return(0);
1817: }
1819: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1820: {
1822: return(0);
1823: }
1825: /*@C
1826: MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
1828: Logically Collective on Mat
1830: Input Parameter:
1831: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1833: Output Parameter:
1834: . lda - the leading dimension
1836: Level: intermediate
1838: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA()
1839: @*/
1840: PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
1841: {
1845: PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1846: return(0);
1847: }
1849: /*@C
1850: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1852: Logically Collective on Mat
1854: Input Parameter:
1855: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1857: Output Parameter:
1858: . array - pointer to the data
1860: Level: intermediate
1862: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1863: @*/
1864: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1865: {
1869: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1870: return(0);
1871: }
1873: /*@C
1874: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1876: Logically Collective on Mat
1878: Input Parameters:
1879: + mat - a MATSEQDENSE or MATMPIDENSE matrix
1880: - array - pointer to the data
1882: Level: intermediate
1884: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1885: @*/
1886: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1887: {
1891: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1892: if (array) *array = NULL;
1893: PetscObjectStateIncrease((PetscObject)A);
1894: return(0);
1895: }
1897: /*@C
1898: MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored
1900: Not Collective
1902: Input Parameter:
1903: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1905: Output Parameter:
1906: . array - pointer to the data
1908: Level: intermediate
1910: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1911: @*/
1912: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1913: {
1917: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1918: return(0);
1919: }
1921: /*@C
1922: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1924: Not Collective
1926: Input Parameters:
1927: + mat - a MATSEQDENSE or MATMPIDENSE matrix
1928: - array - pointer to the data
1930: Level: intermediate
1932: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1933: @*/
1934: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1935: {
1939: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1940: if (array) *array = NULL;
1941: return(0);
1942: }
1944: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1945: {
1946: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1948: PetscInt i,j,nrows,ncols,blda;
1949: const PetscInt *irow,*icol;
1950: PetscScalar *av,*bv,*v = mat->v;
1951: Mat newmat;
1954: ISGetIndices(isrow,&irow);
1955: ISGetIndices(iscol,&icol);
1956: ISGetLocalSize(isrow,&nrows);
1957: ISGetLocalSize(iscol,&ncols);
1959: /* Check submatrixcall */
1960: if (scall == MAT_REUSE_MATRIX) {
1961: PetscInt n_cols,n_rows;
1962: MatGetSize(*B,&n_rows,&n_cols);
1963: if (n_rows != nrows || n_cols != ncols) {
1964: /* resize the result matrix to match number of requested rows/columns */
1965: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1966: }
1967: newmat = *B;
1968: } else {
1969: /* Create and fill new matrix */
1970: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1971: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1972: MatSetType(newmat,((PetscObject)A)->type_name);
1973: MatSeqDenseSetPreallocation(newmat,NULL);
1974: }
1976: /* Now extract the data pointers and do the copy,column at a time */
1977: MatDenseGetArray(newmat,&bv);
1978: MatDenseGetLDA(newmat,&blda);
1979: for (i=0; i<ncols; i++) {
1980: av = v + mat->lda*icol[i];
1981: for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
1982: bv += blda;
1983: }
1984: MatDenseRestoreArray(newmat,&bv);
1986: /* Assemble the matrices so that the correct flags are set */
1987: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1988: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1990: /* Free work space */
1991: ISRestoreIndices(isrow,&irow);
1992: ISRestoreIndices(iscol,&icol);
1993: *B = newmat;
1994: return(0);
1995: }
1997: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1998: {
2000: PetscInt i;
2003: if (scall == MAT_INITIAL_MATRIX) {
2004: PetscCalloc1(n+1,B);
2005: }
2007: for (i=0; i<n; i++) {
2008: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2009: }
2010: return(0);
2011: }
2013: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2014: {
2016: return(0);
2017: }
2019: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2020: {
2022: return(0);
2023: }
2025: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2026: {
2027: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2028: PetscErrorCode ierr;
2029: const PetscScalar *va;
2030: PetscScalar *vb;
2031: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
2034: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2035: if (A->ops->copy != B->ops->copy) {
2036: MatCopy_Basic(A,B,str);
2037: return(0);
2038: }
2039: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2040: MatDenseGetArrayRead(A,&va);
2041: MatDenseGetArray(B,&vb);
2042: if (lda1>m || lda2>m) {
2043: for (j=0; j<n; j++) {
2044: PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2045: }
2046: } else {
2047: PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2048: }
2049: MatDenseRestoreArray(B,&vb);
2050: MatDenseRestoreArrayRead(A,&va);
2051: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2052: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2053: return(0);
2054: }
2056: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2057: {
2061: MatSeqDenseSetPreallocation(A,0);
2062: return(0);
2063: }
2065: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2066: {
2067: PetscInt i,nz = A->rmap->n*A->cmap->n;
2068: PetscScalar *aa;
2072: MatDenseGetArray(A,&aa);
2073: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2074: MatDenseRestoreArray(A,&aa);
2075: return(0);
2076: }
2078: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2079: {
2080: PetscInt i,nz = A->rmap->n*A->cmap->n;
2081: PetscScalar *aa;
2085: MatDenseGetArray(A,&aa);
2086: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2087: MatDenseRestoreArray(A,&aa);
2088: return(0);
2089: }
2091: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2092: {
2093: PetscInt i,nz = A->rmap->n*A->cmap->n;
2094: PetscScalar *aa;
2098: MatDenseGetArray(A,&aa);
2099: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2100: MatDenseRestoreArray(A,&aa);
2101: return(0);
2102: }
2104: /* ----------------------------------------------------------------*/
2105: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2106: {
2110: if (scall == MAT_INITIAL_MATRIX) {
2111: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2112: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2113: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2114: }
2115: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2116: if ((*C)->ops->matmultnumeric) {
2117: (*(*C)->ops->matmultnumeric)(A,B,*C);
2118: } else {
2119: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2120: }
2121: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2122: return(0);
2123: }
2125: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2126: {
2128: PetscInt m=A->rmap->n,n=B->cmap->n;
2129: Mat Cmat;
2130: PetscBool flg;
2133: MatCreate(PETSC_COMM_SELF,&Cmat);
2134: MatSetSizes(Cmat,m,n,m,n);
2135: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2136: MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2137: MatSeqDenseSetPreallocation(Cmat,NULL);
2138: *C = Cmat;
2139: return(0);
2140: }
2142: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2143: {
2144: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2145: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2146: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2147: PetscBLASInt m,n,k;
2148: const PetscScalar *av,*bv;
2149: PetscScalar *cv;
2150: PetscScalar _DOne=1.0,_DZero=0.0;
2151: PetscErrorCode ierr;
2152: PetscBool flg;
2155: /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2156: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2157: if (flg) {
2158: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2159: (*C->ops->matmultnumeric)(A,B,C);
2160: return(0);
2161: }
2162: PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);
2163: if (flg) {
2164: C->ops->matmultnumeric = MatMatMultNumeric_SeqBAIJ_SeqDense;
2165: (*C->ops->matmultnumeric)(A,B,C);
2166: return(0);
2167: }
2168: PetscBLASIntCast(C->rmap->n,&m);
2169: PetscBLASIntCast(C->cmap->n,&n);
2170: PetscBLASIntCast(A->cmap->n,&k);
2171: if (!m || !n || !k) return(0);
2172: MatDenseGetArrayRead(A,&av);
2173: MatDenseGetArrayRead(B,&bv);
2174: MatDenseGetArray(C,&cv);
2175: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2176: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2177: MatDenseRestoreArrayRead(A,&av);
2178: MatDenseRestoreArrayRead(B,&bv);
2179: MatDenseRestoreArray(C,&cv);
2180: return(0);
2181: }
2183: PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2184: {
2188: if (scall == MAT_INITIAL_MATRIX) {
2189: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2190: MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2191: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2192: }
2193: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2194: MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2195: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2196: return(0);
2197: }
2199: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2200: {
2202: PetscInt m=A->rmap->n,n=B->rmap->n;
2203: Mat Cmat;
2204: PetscBool flg;
2207: MatCreate(PETSC_COMM_SELF,&Cmat);
2208: MatSetSizes(Cmat,m,n,m,n);
2209: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2210: MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2211: MatSeqDenseSetPreallocation(Cmat,NULL);
2212: *C = Cmat;
2213: return(0);
2214: }
2216: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2217: {
2218: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2219: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2220: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2221: PetscBLASInt m,n,k;
2222: PetscScalar _DOne=1.0,_DZero=0.0;
2226: PetscBLASIntCast(C->rmap->n,&m);
2227: PetscBLASIntCast(C->cmap->n,&n);
2228: PetscBLASIntCast(A->cmap->n,&k);
2229: if (!m || !n || !k) return(0);
2230: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2231: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2232: return(0);
2233: }
2235: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2236: {
2240: if (scall == MAT_INITIAL_MATRIX) {
2241: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2242: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2243: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2244: }
2245: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2246: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2247: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2248: return(0);
2249: }
2251: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2252: {
2254: PetscInt m=A->cmap->n,n=B->cmap->n;
2255: Mat Cmat;
2256: PetscBool flg;
2259: MatCreate(PETSC_COMM_SELF,&Cmat);
2260: MatSetSizes(Cmat,m,n,m,n);
2261: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2262: MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);
2263: MatSeqDenseSetPreallocation(Cmat,NULL);
2264: *C = Cmat;
2265: return(0);
2266: }
2268: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2269: {
2270: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2271: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2272: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2273: PetscBLASInt m,n,k;
2274: PetscScalar _DOne=1.0,_DZero=0.0;
2278: PetscBLASIntCast(C->rmap->n,&m);
2279: PetscBLASIntCast(C->cmap->n,&n);
2280: PetscBLASIntCast(A->rmap->n,&k);
2281: if (!m || !n || !k) return(0);
2282: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2283: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2284: return(0);
2285: }
2287: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2288: {
2289: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2290: PetscErrorCode ierr;
2291: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2292: PetscScalar *x;
2293: const PetscScalar *aa;
2296: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2297: VecGetArray(v,&x);
2298: VecGetLocalSize(v,&p);
2299: MatDenseGetArrayRead(A,&aa);
2300: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2301: for (i=0; i<m; i++) {
2302: x[i] = aa[i]; if (idx) idx[i] = 0;
2303: for (j=1; j<n; j++) {
2304: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2305: }
2306: }
2307: MatDenseRestoreArrayRead(A,&aa);
2308: VecRestoreArray(v,&x);
2309: return(0);
2310: }
2312: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2313: {
2314: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2315: PetscErrorCode ierr;
2316: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2317: PetscScalar *x;
2318: PetscReal atmp;
2319: const PetscScalar *aa;
2322: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2323: VecGetArray(v,&x);
2324: VecGetLocalSize(v,&p);
2325: MatDenseGetArrayRead(A,&aa);
2326: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2327: for (i=0; i<m; i++) {
2328: x[i] = PetscAbsScalar(aa[i]);
2329: for (j=1; j<n; j++) {
2330: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2331: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2332: }
2333: }
2334: MatDenseRestoreArrayRead(A,&aa);
2335: VecRestoreArray(v,&x);
2336: return(0);
2337: }
2339: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2340: {
2341: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2342: PetscErrorCode ierr;
2343: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2344: PetscScalar *x;
2345: const PetscScalar *aa;
2348: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2349: MatDenseGetArrayRead(A,&aa);
2350: VecGetArray(v,&x);
2351: VecGetLocalSize(v,&p);
2352: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2353: for (i=0; i<m; i++) {
2354: x[i] = aa[i]; if (idx) idx[i] = 0;
2355: for (j=1; j<n; j++) {
2356: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2357: }
2358: }
2359: VecRestoreArray(v,&x);
2360: MatDenseRestoreArrayRead(A,&aa);
2361: return(0);
2362: }
2364: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2365: {
2366: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2367: PetscErrorCode ierr;
2368: PetscScalar *x;
2369: const PetscScalar *aa;
2372: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2373: MatDenseGetArrayRead(A,&aa);
2374: VecGetArray(v,&x);
2375: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2376: VecRestoreArray(v,&x);
2377: MatDenseRestoreArrayRead(A,&aa);
2378: return(0);
2379: }
2381: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2382: {
2383: PetscErrorCode ierr;
2384: PetscInt i,j,m,n;
2385: const PetscScalar *a;
2388: MatGetSize(A,&m,&n);
2389: PetscArrayzero(norms,n);
2390: MatDenseGetArrayRead(A,&a);
2391: if (type == NORM_2) {
2392: for (i=0; i<n; i++) {
2393: for (j=0; j<m; j++) {
2394: norms[i] += PetscAbsScalar(a[j]*a[j]);
2395: }
2396: a += m;
2397: }
2398: } else if (type == NORM_1) {
2399: for (i=0; i<n; i++) {
2400: for (j=0; j<m; j++) {
2401: norms[i] += PetscAbsScalar(a[j]);
2402: }
2403: a += m;
2404: }
2405: } else if (type == NORM_INFINITY) {
2406: for (i=0; i<n; i++) {
2407: for (j=0; j<m; j++) {
2408: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2409: }
2410: a += m;
2411: }
2412: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2413: MatDenseRestoreArrayRead(A,&a);
2414: if (type == NORM_2) {
2415: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2416: }
2417: return(0);
2418: }
2420: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2421: {
2423: PetscScalar *a;
2424: PetscInt m,n,i;
2427: MatGetSize(x,&m,&n);
2428: MatDenseGetArray(x,&a);
2429: for (i=0; i<m*n; i++) {
2430: PetscRandomGetValue(rctx,a+i);
2431: }
2432: MatDenseRestoreArray(x,&a);
2433: return(0);
2434: }
2436: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2437: {
2439: *missing = PETSC_FALSE;
2440: return(0);
2441: }
2443: /* vals is not const */
2444: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2445: {
2447: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2448: PetscScalar *v;
2451: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2452: MatDenseGetArray(A,&v);
2453: *vals = v+col*a->lda;
2454: MatDenseRestoreArray(A,&v);
2455: return(0);
2456: }
2458: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2459: {
2461: *vals = 0; /* user cannot accidently use the array later */
2462: return(0);
2463: }
2465: /* -------------------------------------------------------------------*/
2466: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2467: MatGetRow_SeqDense,
2468: MatRestoreRow_SeqDense,
2469: MatMult_SeqDense,
2470: /* 4*/ MatMultAdd_SeqDense,
2471: MatMultTranspose_SeqDense,
2472: MatMultTransposeAdd_SeqDense,
2473: 0,
2474: 0,
2475: 0,
2476: /* 10*/ 0,
2477: MatLUFactor_SeqDense,
2478: MatCholeskyFactor_SeqDense,
2479: MatSOR_SeqDense,
2480: MatTranspose_SeqDense,
2481: /* 15*/ MatGetInfo_SeqDense,
2482: MatEqual_SeqDense,
2483: MatGetDiagonal_SeqDense,
2484: MatDiagonalScale_SeqDense,
2485: MatNorm_SeqDense,
2486: /* 20*/ MatAssemblyBegin_SeqDense,
2487: MatAssemblyEnd_SeqDense,
2488: MatSetOption_SeqDense,
2489: MatZeroEntries_SeqDense,
2490: /* 24*/ MatZeroRows_SeqDense,
2491: 0,
2492: 0,
2493: 0,
2494: 0,
2495: /* 29*/ MatSetUp_SeqDense,
2496: 0,
2497: 0,
2498: 0,
2499: 0,
2500: /* 34*/ MatDuplicate_SeqDense,
2501: 0,
2502: 0,
2503: 0,
2504: 0,
2505: /* 39*/ MatAXPY_SeqDense,
2506: MatCreateSubMatrices_SeqDense,
2507: 0,
2508: MatGetValues_SeqDense,
2509: MatCopy_SeqDense,
2510: /* 44*/ MatGetRowMax_SeqDense,
2511: MatScale_SeqDense,
2512: MatShift_Basic,
2513: 0,
2514: MatZeroRowsColumns_SeqDense,
2515: /* 49*/ MatSetRandom_SeqDense,
2516: 0,
2517: 0,
2518: 0,
2519: 0,
2520: /* 54*/ 0,
2521: 0,
2522: 0,
2523: 0,
2524: 0,
2525: /* 59*/ 0,
2526: MatDestroy_SeqDense,
2527: MatView_SeqDense,
2528: 0,
2529: 0,
2530: /* 64*/ 0,
2531: 0,
2532: 0,
2533: 0,
2534: 0,
2535: /* 69*/ MatGetRowMaxAbs_SeqDense,
2536: 0,
2537: 0,
2538: 0,
2539: 0,
2540: /* 74*/ 0,
2541: 0,
2542: 0,
2543: 0,
2544: 0,
2545: /* 79*/ 0,
2546: 0,
2547: 0,
2548: 0,
2549: /* 83*/ MatLoad_SeqDense,
2550: 0,
2551: MatIsHermitian_SeqDense,
2552: 0,
2553: 0,
2554: 0,
2555: /* 89*/ MatMatMult_SeqDense_SeqDense,
2556: MatMatMultSymbolic_SeqDense_SeqDense,
2557: MatMatMultNumeric_SeqDense_SeqDense,
2558: MatPtAP_SeqDense_SeqDense,
2559: MatPtAPSymbolic_SeqDense_SeqDense,
2560: /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2561: MatMatTransposeMult_SeqDense_SeqDense,
2562: MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2563: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2564: 0,
2565: /* 99*/ 0,
2566: 0,
2567: 0,
2568: MatConjugate_SeqDense,
2569: 0,
2570: /*104*/ 0,
2571: MatRealPart_SeqDense,
2572: MatImaginaryPart_SeqDense,
2573: 0,
2574: 0,
2575: /*109*/ 0,
2576: 0,
2577: MatGetRowMin_SeqDense,
2578: MatGetColumnVector_SeqDense,
2579: MatMissingDiagonal_SeqDense,
2580: /*114*/ 0,
2581: 0,
2582: 0,
2583: 0,
2584: 0,
2585: /*119*/ 0,
2586: 0,
2587: 0,
2588: 0,
2589: 0,
2590: /*124*/ 0,
2591: MatGetColumnNorms_SeqDense,
2592: 0,
2593: 0,
2594: 0,
2595: /*129*/ 0,
2596: MatTransposeMatMult_SeqDense_SeqDense,
2597: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2598: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2599: 0,
2600: /*134*/ 0,
2601: 0,
2602: 0,
2603: 0,
2604: 0,
2605: /*139*/ 0,
2606: 0,
2607: 0,
2608: 0,
2609: 0,
2610: /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2611: };
2613: /*@C
2614: MatCreateSeqDense - Creates a sequential dense matrix that
2615: is stored in column major order (the usual Fortran 77 manner). Many
2616: of the matrix operations use the BLAS and LAPACK routines.
2618: Collective
2620: Input Parameters:
2621: + comm - MPI communicator, set to PETSC_COMM_SELF
2622: . m - number of rows
2623: . n - number of columns
2624: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2625: to control all matrix memory allocation.
2627: Output Parameter:
2628: . A - the matrix
2630: Notes:
2631: The data input variable is intended primarily for Fortran programmers
2632: who wish to allocate their own matrix memory space. Most users should
2633: set data=NULL.
2635: Level: intermediate
2637: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2638: @*/
2639: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2640: {
2644: MatCreate(comm,A);
2645: MatSetSizes(*A,m,n,m,n);
2646: MatSetType(*A,MATSEQDENSE);
2647: MatSeqDenseSetPreallocation(*A,data);
2648: return(0);
2649: }
2651: /*@C
2652: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2654: Collective
2656: Input Parameters:
2657: + B - the matrix
2658: - data - the array (or NULL)
2660: Notes:
2661: The data input variable is intended primarily for Fortran programmers
2662: who wish to allocate their own matrix memory space. Most users should
2663: need not call this routine.
2665: Level: intermediate
2667: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2669: @*/
2670: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2671: {
2675: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2676: return(0);
2677: }
2679: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2680: {
2681: Mat_SeqDense *b;
2685: B->preallocated = PETSC_TRUE;
2687: PetscLayoutSetUp(B->rmap);
2688: PetscLayoutSetUp(B->cmap);
2690: b = (Mat_SeqDense*)B->data;
2691: b->Mmax = B->rmap->n;
2692: b->Nmax = B->cmap->n;
2693: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2695: PetscIntMultError(b->lda,b->Nmax,NULL);
2696: if (!data) { /* petsc-allocated storage */
2697: if (!b->user_alloc) { PetscFree(b->v); }
2698: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2699: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2701: b->user_alloc = PETSC_FALSE;
2702: } else { /* user-allocated storage */
2703: if (!b->user_alloc) { PetscFree(b->v); }
2704: b->v = data;
2705: b->user_alloc = PETSC_TRUE;
2706: }
2707: B->assembled = PETSC_TRUE;
2708: return(0);
2709: }
2711: #if defined(PETSC_HAVE_ELEMENTAL)
2712: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2713: {
2714: Mat mat_elemental;
2715: PetscErrorCode ierr;
2716: const PetscScalar *array;
2717: PetscScalar *v_colwise;
2718: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2721: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2722: MatDenseGetArrayRead(A,&array);
2723: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2724: k = 0;
2725: for (j=0; j<N; j++) {
2726: cols[j] = j;
2727: for (i=0; i<M; i++) {
2728: v_colwise[j*M+i] = array[k++];
2729: }
2730: }
2731: for (i=0; i<M; i++) {
2732: rows[i] = i;
2733: }
2734: MatDenseRestoreArrayRead(A,&array);
2736: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2737: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2738: MatSetType(mat_elemental,MATELEMENTAL);
2739: MatSetUp(mat_elemental);
2741: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2742: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2743: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2744: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2745: PetscFree3(v_colwise,rows,cols);
2747: if (reuse == MAT_INPLACE_MATRIX) {
2748: MatHeaderReplace(A,&mat_elemental);
2749: } else {
2750: *newmat = mat_elemental;
2751: }
2752: return(0);
2753: }
2754: #endif
2756: /*@C
2757: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2759: Input parameter:
2760: + A - the matrix
2761: - lda - the leading dimension
2763: Notes:
2764: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2765: it asserts that the preallocation has a leading dimension (the LDA parameter
2766: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2768: Level: intermediate
2770: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2772: @*/
2773: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2774: {
2775: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2778: if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n);
2779: b->lda = lda;
2780: b->changelda = PETSC_FALSE;
2781: b->Mmax = PetscMax(b->Mmax,lda);
2782: return(0);
2783: }
2785: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2786: {
2788: PetscMPIInt size;
2791: MPI_Comm_size(comm,&size);
2792: if (size == 1) {
2793: if (scall == MAT_INITIAL_MATRIX) {
2794: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2795: } else {
2796: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2797: }
2798: } else {
2799: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2800: }
2801: return(0);
2802: }
2804: /*MC
2805: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2807: Options Database Keys:
2808: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2810: Level: beginner
2812: .seealso: MatCreateSeqDense()
2814: M*/
2815: PetscErrorCode MatCreate_SeqDense(Mat B)
2816: {
2817: Mat_SeqDense *b;
2819: PetscMPIInt size;
2822: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2823: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2825: PetscNewLog(B,&b);
2826: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2827: B->data = (void*)b;
2829: b->roworiented = PETSC_TRUE;
2831: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
2832: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2833: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2834: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2835: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2836: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2837: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2838: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2839: #if defined(PETSC_HAVE_ELEMENTAL)
2840: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2841: #endif
2842: #if defined(PETSC_HAVE_CUDA)
2843: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
2844: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijcusparse_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2845: #endif
2846: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2847: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2848: #if defined(PETSC_HAVE_VIENNACL)
2849: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijviennacl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2850: #endif
2851: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2852: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2853: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqbaij_seqdense_C",MatMatMult_SeqBAIJ_SeqDense);
2854: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);
2855: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);
2856: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2857: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2858: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2859: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2860: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2861: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2862: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2863: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2864: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);
2865: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2866: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2867: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2868: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);
2870: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2871: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2872: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2873: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2874: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2875: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2876: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2877: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2878: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2880: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2881: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2882: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2883: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2884: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2885: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2886: return(0);
2887: }
2889: /*@C
2890: 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.
2892: Not Collective
2894: Input Parameter:
2895: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2896: - col - column index
2898: Output Parameter:
2899: . vals - pointer to the data
2901: Level: intermediate
2903: .seealso: MatDenseRestoreColumn()
2904: @*/
2905: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2906: {
2910: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2911: return(0);
2912: }
2914: /*@C
2915: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
2917: Not Collective
2919: Input Parameter:
2920: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2922: Output Parameter:
2923: . vals - pointer to the data
2925: Level: intermediate
2927: .seealso: MatDenseGetColumn()
2928: @*/
2929: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2930: {
2934: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2935: return(0);
2936: }