Actual source code: dense.c
petsc-3.7.5 2017-01-01
2: /*
3: Defines the basic matrix operations for sequential dense.
4: */
6: #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/
7: #include <petscblaslapack.h>
9: #include <../src/mat/impls/aij/seq/aij.h>
13: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
14: {
15: PetscErrorCode ierr;
16: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
17: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
18: PetscScalar *slot,*bb;
19: const PetscScalar *xx;
22: #if defined(PETSC_USE_DEBUG)
23: for (i=0; i<N; i++) {
24: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
25: 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);
26: 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);
27: }
28: #endif
30: /* fix right hand side if needed */
31: if (x && b) {
32: VecGetArrayRead(x,&xx);
33: VecGetArray(b,&bb);
34: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
35: VecRestoreArrayRead(x,&xx);
36: VecRestoreArray(b,&bb);
37: }
39: for (i=0; i<N; i++) {
40: slot = l->v + rows[i]*m;
41: PetscMemzero(slot,r*sizeof(PetscScalar));
42: }
43: for (i=0; i<N; i++) {
44: slot = l->v + rows[i];
45: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
46: }
47: if (diag != 0.0) {
48: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
49: for (i=0; i<N; i++) {
50: slot = l->v + (m+1)*rows[i];
51: *slot = diag;
52: }
53: }
54: return(0);
55: }
59: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
60: {
61: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
65: MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
66: MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
67: return(0);
68: }
72: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
73: {
74: Mat_SeqDense *c;
78: MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
79: c = (Mat_SeqDense*)((*C)->data);
80: MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
81: return(0);
82: }
86: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
87: {
91: if (reuse == MAT_INITIAL_MATRIX) {
92: MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
93: }
94: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
95: (*(*C)->ops->ptapnumeric)(A,P,*C);
96: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
97: return(0);
98: }
102: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
103: {
104: Mat B;
105: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
106: Mat_SeqDense *b;
108: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
109: MatScalar *av=a->a;
112: MatCreate(PetscObjectComm((PetscObject)A),&B);
113: MatSetSizes(B,m,n,m,n);
114: MatSetType(B,MATSEQDENSE);
115: MatSeqDenseSetPreallocation(B,NULL);
117: b = (Mat_SeqDense*)(B->data);
118: for (i=0; i<m; i++) {
119: PetscInt j;
120: for (j=0;j<ai[1]-ai[0];j++) {
121: b->v[*aj*m+i] = *av;
122: aj++;
123: av++;
124: }
125: ai++;
126: }
127: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
130: if (reuse == MAT_INPLACE_MATRIX) {
131: MatHeaderReplace(A,&B);
132: } else {
133: *newmat = B;
134: }
135: return(0);
136: }
140: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
141: {
142: Mat B;
143: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
145: PetscInt i, j;
146: PetscInt *rows, *nnz;
147: MatScalar *aa = a->v, *vals;
150: MatCreate(PetscObjectComm((PetscObject)A),&B);
151: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
152: MatSetType(B,MATSEQAIJ);
153: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
154: for (j=0; j<A->cmap->n; j++) {
155: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
156: aa += a->lda;
157: }
158: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
159: aa = a->v;
160: for (j=0; j<A->cmap->n; j++) {
161: PetscInt numRows = 0;
162: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
163: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
164: aa += a->lda;
165: }
166: PetscFree3(rows,nnz,vals);
167: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
170: if (reuse == MAT_INPLACE_MATRIX) {
171: MatHeaderReplace(A,&B);
172: } else {
173: *newmat = B;
174: }
175: return(0);
176: }
180: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
181: {
182: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
183: PetscScalar oalpha = alpha;
184: PetscInt j;
185: PetscBLASInt N,m,ldax,lday,one = 1;
189: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
190: PetscBLASIntCast(X->rmap->n,&m);
191: PetscBLASIntCast(x->lda,&ldax);
192: PetscBLASIntCast(y->lda,&lday);
193: if (ldax>m || lday>m) {
194: for (j=0; j<X->cmap->n; j++) {
195: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
196: }
197: } else {
198: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
199: }
200: PetscObjectStateIncrease((PetscObject)Y);
201: PetscLogFlops(PetscMax(2*N-1,0));
202: return(0);
203: }
207: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
208: {
209: PetscInt N = A->rmap->n*A->cmap->n;
212: info->block_size = 1.0;
213: info->nz_allocated = (double)N;
214: info->nz_used = (double)N;
215: info->nz_unneeded = (double)0;
216: info->assemblies = (double)A->num_ass;
217: info->mallocs = 0;
218: info->memory = ((PetscObject)A)->mem;
219: info->fill_ratio_given = 0;
220: info->fill_ratio_needed = 0;
221: info->factor_mallocs = 0;
222: return(0);
223: }
227: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
228: {
229: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
230: PetscScalar oalpha = alpha;
232: PetscBLASInt one = 1,j,nz,lda;
235: PetscBLASIntCast(a->lda,&lda);
236: if (lda>A->rmap->n) {
237: PetscBLASIntCast(A->rmap->n,&nz);
238: for (j=0; j<A->cmap->n; j++) {
239: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
240: }
241: } else {
242: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
243: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
244: }
245: PetscLogFlops(nz);
246: return(0);
247: }
251: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
252: {
253: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
254: PetscInt i,j,m = A->rmap->n,N;
255: PetscScalar *v = a->v;
258: *fl = PETSC_FALSE;
259: if (A->rmap->n != A->cmap->n) return(0);
260: N = a->lda;
262: for (i=0; i<m; i++) {
263: for (j=i+1; j<m; j++) {
264: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
265: }
266: }
267: *fl = PETSC_TRUE;
268: return(0);
269: }
273: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
274: {
275: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
277: PetscInt lda = (PetscInt)mat->lda,j,m;
280: PetscLayoutReference(A->rmap,&newi->rmap);
281: PetscLayoutReference(A->cmap,&newi->cmap);
282: MatSeqDenseSetPreallocation(newi,NULL);
283: if (cpvalues == MAT_COPY_VALUES) {
284: l = (Mat_SeqDense*)newi->data;
285: if (lda>A->rmap->n) {
286: m = A->rmap->n;
287: for (j=0; j<A->cmap->n; j++) {
288: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
289: }
290: } else {
291: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
292: }
293: }
294: newi->assembled = PETSC_TRUE;
295: return(0);
296: }
300: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
301: {
305: MatCreate(PetscObjectComm((PetscObject)A),newmat);
306: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
307: MatSetType(*newmat,((PetscObject)A)->type_name);
308: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
309: return(0);
310: }
313: static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
317: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
318: {
319: MatFactorInfo info;
323: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
324: MatLUFactor_SeqDense(fact,0,0,&info);
325: return(0);
326: }
330: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
331: {
332: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
333: PetscErrorCode ierr;
334: const PetscScalar *x;
335: PetscScalar *y;
336: PetscBLASInt one = 1,info,m;
339: PetscBLASIntCast(A->rmap->n,&m);
340: VecGetArrayRead(xx,&x);
341: VecGetArray(yy,&y);
342: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
343: if (A->factortype == MAT_FACTOR_LU) {
344: #if defined(PETSC_MISSING_LAPACK_GETRS)
345: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
346: #else
347: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
348: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
349: #endif
350: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
351: #if defined(PETSC_MISSING_LAPACK_POTRS)
352: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
353: #else
354: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
355: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
356: #endif
357: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
358: VecRestoreArrayRead(xx,&x);
359: VecRestoreArray(yy,&y);
360: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
361: return(0);
362: }
366: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
367: {
368: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
370: PetscScalar *b,*x;
371: PetscInt n;
372: PetscBLASInt nrhs,info,m;
373: PetscBool flg;
376: PetscBLASIntCast(A->rmap->n,&m);
377: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
378: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
379: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
380: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
382: MatGetSize(B,NULL,&n);
383: PetscBLASIntCast(n,&nrhs);
384: MatDenseGetArray(B,&b);
385: MatDenseGetArray(X,&x);
387: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
389: if (A->factortype == MAT_FACTOR_LU) {
390: #if defined(PETSC_MISSING_LAPACK_GETRS)
391: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
392: #else
393: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
394: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
395: #endif
396: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
397: #if defined(PETSC_MISSING_LAPACK_POTRS)
398: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
399: #else
400: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
401: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
402: #endif
403: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
405: MatDenseRestoreArray(B,&b);
406: MatDenseRestoreArray(X,&x);
407: PetscLogFlops(nrhs*(2.0*m*m - m));
408: return(0);
409: }
413: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
414: {
415: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
416: PetscErrorCode ierr;
417: const PetscScalar *x;
418: PetscScalar *y;
419: PetscBLASInt one = 1,info,m;
422: PetscBLASIntCast(A->rmap->n,&m);
423: VecGetArrayRead(xx,&x);
424: VecGetArray(yy,&y);
425: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
426: /* assume if pivots exist then use LU; else Cholesky */
427: if (mat->pivots) {
428: #if defined(PETSC_MISSING_LAPACK_GETRS)
429: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
430: #else
431: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
432: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
433: #endif
434: } else {
435: #if defined(PETSC_MISSING_LAPACK_POTRS)
436: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
437: #else
438: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
439: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
440: #endif
441: }
442: VecRestoreArrayRead(xx,&x);
443: VecRestoreArray(yy,&y);
444: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
445: return(0);
446: }
450: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
451: {
452: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
453: PetscErrorCode ierr;
454: const PetscScalar *x;
455: PetscScalar *y,sone = 1.0;
456: Vec tmp = 0;
457: PetscBLASInt one = 1,info,m;
460: PetscBLASIntCast(A->rmap->n,&m);
461: VecGetArrayRead(xx,&x);
462: VecGetArray(yy,&y);
463: if (!A->rmap->n || !A->cmap->n) return(0);
464: if (yy == zz) {
465: VecDuplicate(yy,&tmp);
466: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
467: VecCopy(yy,tmp);
468: }
469: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
470: /* assume if pivots exist then use LU; else Cholesky */
471: if (mat->pivots) {
472: #if defined(PETSC_MISSING_LAPACK_GETRS)
473: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
474: #else
475: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
476: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
477: #endif
478: } else {
479: #if defined(PETSC_MISSING_LAPACK_POTRS)
480: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
481: #else
482: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
483: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
484: #endif
485: }
486: if (tmp) {
487: VecAXPY(yy,sone,tmp);
488: VecDestroy(&tmp);
489: } else {
490: VecAXPY(yy,sone,zz);
491: }
492: VecRestoreArrayRead(xx,&x);
493: VecRestoreArray(yy,&y);
494: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
495: return(0);
496: }
500: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
501: {
502: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
503: PetscErrorCode ierr;
504: const PetscScalar *x;
505: PetscScalar *y,sone = 1.0;
506: Vec tmp = NULL;
507: PetscBLASInt one = 1,info,m;
510: PetscBLASIntCast(A->rmap->n,&m);
511: if (!A->rmap->n || !A->cmap->n) return(0);
512: VecGetArrayRead(xx,&x);
513: VecGetArray(yy,&y);
514: if (yy == zz) {
515: VecDuplicate(yy,&tmp);
516: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
517: VecCopy(yy,tmp);
518: }
519: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
520: /* assume if pivots exist then use LU; else Cholesky */
521: if (mat->pivots) {
522: #if defined(PETSC_MISSING_LAPACK_GETRS)
523: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
524: #else
525: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
526: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
527: #endif
528: } else {
529: #if defined(PETSC_MISSING_LAPACK_POTRS)
530: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
531: #else
532: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
533: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
534: #endif
535: }
536: if (tmp) {
537: VecAXPY(yy,sone,tmp);
538: VecDestroy(&tmp);
539: } else {
540: VecAXPY(yy,sone,zz);
541: }
542: VecRestoreArrayRead(xx,&x);
543: VecRestoreArray(yy,&y);
544: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
545: return(0);
546: }
548: /* ---------------------------------------------------------------*/
549: /* COMMENT: I have chosen to hide row permutation in the pivots,
550: rather than put it in the Mat->row slot.*/
553: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
554: {
555: #if defined(PETSC_MISSING_LAPACK_GETRF)
557: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
558: #else
559: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
561: PetscBLASInt n,m,info;
564: PetscBLASIntCast(A->cmap->n,&n);
565: PetscBLASIntCast(A->rmap->n,&m);
566: if (!mat->pivots) {
567: PetscMalloc1(A->rmap->n+1,&mat->pivots);
568: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
569: }
570: if (!A->rmap->n || !A->cmap->n) return(0);
571: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
572: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
573: PetscFPTrapPop();
575: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
576: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
577: A->ops->solve = MatSolve_SeqDense;
578: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
579: A->ops->solveadd = MatSolveAdd_SeqDense;
580: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
581: A->factortype = MAT_FACTOR_LU;
583: PetscFree(A->solvertype);
584: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
586: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
587: #endif
588: return(0);
589: }
593: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
594: {
595: #if defined(PETSC_MISSING_LAPACK_POTRF)
597: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
598: #else
599: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
601: PetscBLASInt info,n;
604: PetscBLASIntCast(A->cmap->n,&n);
605: PetscFree(mat->pivots);
607: if (!A->rmap->n || !A->cmap->n) return(0);
608: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
609: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
610: A->ops->solve = MatSolve_SeqDense;
611: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
612: A->ops->solveadd = MatSolveAdd_SeqDense;
613: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
614: A->factortype = MAT_FACTOR_CHOLESKY;
616: PetscFree(A->solvertype);
617: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
619: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
620: #endif
621: return(0);
622: }
627: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
628: {
630: MatFactorInfo info;
633: info.fill = 1.0;
635: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
636: MatCholeskyFactor_SeqDense(fact,0,&info);
637: return(0);
638: }
642: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
643: {
645: fact->assembled = PETSC_TRUE;
646: fact->preallocated = PETSC_TRUE;
647: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
648: return(0);
649: }
653: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
654: {
656: fact->preallocated = PETSC_TRUE;
657: fact->assembled = PETSC_TRUE;
658: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
659: return(0);
660: }
664: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
665: {
669: MatCreate(PetscObjectComm((PetscObject)A),fact);
670: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
671: MatSetType(*fact,((PetscObject)A)->type_name);
672: if (ftype == MAT_FACTOR_LU) {
673: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
674: } else {
675: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
676: }
677: (*fact)->factortype = ftype;
679: PetscFree((*fact)->solvertype);
680: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
681: return(0);
682: }
684: /* ------------------------------------------------------------------*/
687: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
688: {
689: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
690: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
691: const PetscScalar *b;
692: PetscErrorCode ierr;
693: PetscInt m = A->rmap->n,i;
694: PetscBLASInt o = 1,bm;
697: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
698: PetscBLASIntCast(m,&bm);
699: if (flag & SOR_ZERO_INITIAL_GUESS) {
700: /* this is a hack fix, should have another version without the second BLASdot */
701: VecSet(xx,zero);
702: }
703: VecGetArray(xx,&x);
704: VecGetArrayRead(bb,&b);
705: its = its*lits;
706: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
707: while (its--) {
708: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
709: for (i=0; i<m; i++) {
710: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
711: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
712: }
713: }
714: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
715: for (i=m-1; i>=0; i--) {
716: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
717: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
718: }
719: }
720: }
721: VecRestoreArrayRead(bb,&b);
722: VecRestoreArray(xx,&x);
723: return(0);
724: }
726: /* -----------------------------------------------------------------*/
729: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
730: {
731: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
732: const PetscScalar *v = mat->v,*x;
733: PetscScalar *y;
734: PetscErrorCode ierr;
735: PetscBLASInt m, n,_One=1;
736: PetscScalar _DOne=1.0,_DZero=0.0;
739: PetscBLASIntCast(A->rmap->n,&m);
740: PetscBLASIntCast(A->cmap->n,&n);
741: if (!A->rmap->n || !A->cmap->n) return(0);
742: VecGetArrayRead(xx,&x);
743: VecGetArray(yy,&y);
744: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
745: VecRestoreArrayRead(xx,&x);
746: VecRestoreArray(yy,&y);
747: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
748: return(0);
749: }
753: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
754: {
755: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
756: PetscScalar *y,_DOne=1.0,_DZero=0.0;
757: PetscErrorCode ierr;
758: PetscBLASInt m, n, _One=1;
759: const PetscScalar *v = mat->v,*x;
762: PetscBLASIntCast(A->rmap->n,&m);
763: PetscBLASIntCast(A->cmap->n,&n);
764: if (!A->rmap->n || !A->cmap->n) return(0);
765: VecGetArrayRead(xx,&x);
766: VecGetArray(yy,&y);
767: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
768: VecRestoreArrayRead(xx,&x);
769: VecRestoreArray(yy,&y);
770: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
771: return(0);
772: }
776: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
777: {
778: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
779: const PetscScalar *v = mat->v,*x;
780: PetscScalar *y,_DOne=1.0;
781: PetscErrorCode ierr;
782: PetscBLASInt m, n, _One=1;
785: PetscBLASIntCast(A->rmap->n,&m);
786: PetscBLASIntCast(A->cmap->n,&n);
787: if (!A->rmap->n || !A->cmap->n) return(0);
788: if (zz != yy) {VecCopy(zz,yy);}
789: VecGetArrayRead(xx,&x);
790: VecGetArray(yy,&y);
791: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
792: VecRestoreArrayRead(xx,&x);
793: VecRestoreArray(yy,&y);
794: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
795: return(0);
796: }
800: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
801: {
802: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
803: const PetscScalar *v = mat->v,*x;
804: PetscScalar *y;
805: PetscErrorCode ierr;
806: PetscBLASInt m, n, _One=1;
807: PetscScalar _DOne=1.0;
810: PetscBLASIntCast(A->rmap->n,&m);
811: PetscBLASIntCast(A->cmap->n,&n);
812: if (!A->rmap->n || !A->cmap->n) return(0);
813: if (zz != yy) {VecCopy(zz,yy);}
814: VecGetArrayRead(xx,&x);
815: VecGetArray(yy,&y);
816: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
817: VecRestoreArrayRead(xx,&x);
818: VecRestoreArray(yy,&y);
819: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
820: return(0);
821: }
823: /* -----------------------------------------------------------------*/
826: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
827: {
828: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
829: PetscScalar *v;
831: PetscInt i;
834: *ncols = A->cmap->n;
835: if (cols) {
836: PetscMalloc1(A->cmap->n+1,cols);
837: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
838: }
839: if (vals) {
840: PetscMalloc1(A->cmap->n+1,vals);
841: v = mat->v + row;
842: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
843: }
844: return(0);
845: }
849: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
850: {
854: if (cols) {PetscFree(*cols);}
855: if (vals) {PetscFree(*vals); }
856: return(0);
857: }
858: /* ----------------------------------------------------------------*/
861: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
862: {
863: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
864: PetscInt i,j,idx=0;
867: if (!mat->roworiented) {
868: if (addv == INSERT_VALUES) {
869: for (j=0; j<n; j++) {
870: if (indexn[j] < 0) {idx += m; continue;}
871: #if defined(PETSC_USE_DEBUG)
872: 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);
873: #endif
874: for (i=0; i<m; i++) {
875: if (indexm[i] < 0) {idx++; continue;}
876: #if defined(PETSC_USE_DEBUG)
877: 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);
878: #endif
879: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
880: }
881: }
882: } else {
883: for (j=0; j<n; j++) {
884: if (indexn[j] < 0) {idx += m; continue;}
885: #if defined(PETSC_USE_DEBUG)
886: 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);
887: #endif
888: for (i=0; i<m; i++) {
889: if (indexm[i] < 0) {idx++; continue;}
890: #if defined(PETSC_USE_DEBUG)
891: 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);
892: #endif
893: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
894: }
895: }
896: }
897: } else {
898: if (addv == INSERT_VALUES) {
899: for (i=0; i<m; i++) {
900: if (indexm[i] < 0) { idx += n; continue;}
901: #if defined(PETSC_USE_DEBUG)
902: 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);
903: #endif
904: for (j=0; j<n; j++) {
905: if (indexn[j] < 0) { idx++; continue;}
906: #if defined(PETSC_USE_DEBUG)
907: 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);
908: #endif
909: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
910: }
911: }
912: } else {
913: for (i=0; i<m; i++) {
914: if (indexm[i] < 0) { idx += n; continue;}
915: #if defined(PETSC_USE_DEBUG)
916: 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);
917: #endif
918: for (j=0; j<n; j++) {
919: if (indexn[j] < 0) { idx++; continue;}
920: #if defined(PETSC_USE_DEBUG)
921: 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);
922: #endif
923: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
924: }
925: }
926: }
927: }
928: return(0);
929: }
933: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
934: {
935: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
936: PetscInt i,j;
939: /* row-oriented output */
940: for (i=0; i<m; i++) {
941: if (indexm[i] < 0) {v += n;continue;}
942: 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);
943: for (j=0; j<n; j++) {
944: if (indexn[j] < 0) {v++; continue;}
945: 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);
946: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
947: }
948: }
949: return(0);
950: }
952: /* -----------------------------------------------------------------*/
956: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
957: {
958: Mat_SeqDense *a;
960: PetscInt *scols,i,j,nz,header[4];
961: int fd;
962: PetscMPIInt size;
963: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
964: PetscScalar *vals,*svals,*v,*w;
965: MPI_Comm comm;
968: /* force binary viewer to load .info file if it has not yet done so */
969: PetscViewerSetUp(viewer);
970: PetscObjectGetComm((PetscObject)viewer,&comm);
971: MPI_Comm_size(comm,&size);
972: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
973: PetscViewerBinaryGetDescriptor(viewer,&fd);
974: PetscBinaryRead(fd,header,4,PETSC_INT);
975: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
976: M = header[1]; N = header[2]; nz = header[3];
978: /* set global size if not set already*/
979: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
980: MatSetSizes(newmat,M,N,M,N);
981: } else {
982: /* if sizes and type are already set, check if the vector global sizes are correct */
983: MatGetSize(newmat,&grows,&gcols);
984: 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);
985: }
986: a = (Mat_SeqDense*)newmat->data;
987: if (!a->user_alloc) {
988: MatSeqDenseSetPreallocation(newmat,NULL);
989: }
991: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
992: a = (Mat_SeqDense*)newmat->data;
993: v = a->v;
994: /* Allocate some temp space to read in the values and then flip them
995: from row major to column major */
996: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
997: /* read in nonzero values */
998: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
999: /* now flip the values and store them in the matrix*/
1000: for (j=0; j<N; j++) {
1001: for (i=0; i<M; i++) {
1002: *v++ =w[i*N+j];
1003: }
1004: }
1005: PetscFree(w);
1006: } else {
1007: /* read row lengths */
1008: PetscMalloc1(M+1,&rowlengths);
1009: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1011: a = (Mat_SeqDense*)newmat->data;
1012: v = a->v;
1014: /* read column indices and nonzeros */
1015: PetscMalloc1(nz+1,&scols);
1016: cols = scols;
1017: PetscBinaryRead(fd,cols,nz,PETSC_INT);
1018: PetscMalloc1(nz+1,&svals);
1019: vals = svals;
1020: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
1022: /* insert into matrix */
1023: for (i=0; i<M; i++) {
1024: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1025: svals += rowlengths[i]; scols += rowlengths[i];
1026: }
1027: PetscFree(vals);
1028: PetscFree(cols);
1029: PetscFree(rowlengths);
1030: }
1031: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1032: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1033: return(0);
1034: }
1038: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1039: {
1040: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1041: PetscErrorCode ierr;
1042: PetscInt i,j;
1043: const char *name;
1044: PetscScalar *v;
1045: PetscViewerFormat format;
1046: #if defined(PETSC_USE_COMPLEX)
1047: PetscBool allreal = PETSC_TRUE;
1048: #endif
1051: PetscViewerGetFormat(viewer,&format);
1052: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1053: return(0); /* do nothing for now */
1054: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1055: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1056: for (i=0; i<A->rmap->n; i++) {
1057: v = a->v + i;
1058: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1059: for (j=0; j<A->cmap->n; j++) {
1060: #if defined(PETSC_USE_COMPLEX)
1061: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1062: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1063: } else if (PetscRealPart(*v)) {
1064: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1065: }
1066: #else
1067: if (*v) {
1068: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1069: }
1070: #endif
1071: v += a->lda;
1072: }
1073: PetscViewerASCIIPrintf(viewer,"\n");
1074: }
1075: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1076: } else {
1077: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1078: #if defined(PETSC_USE_COMPLEX)
1079: /* determine if matrix has all real values */
1080: v = a->v;
1081: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1082: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1083: }
1084: #endif
1085: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1086: PetscObjectGetName((PetscObject)A,&name);
1087: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1088: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1089: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1090: }
1092: for (i=0; i<A->rmap->n; i++) {
1093: v = a->v + i;
1094: for (j=0; j<A->cmap->n; j++) {
1095: #if defined(PETSC_USE_COMPLEX)
1096: if (allreal) {
1097: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1098: } else {
1099: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1100: }
1101: #else
1102: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1103: #endif
1104: v += a->lda;
1105: }
1106: PetscViewerASCIIPrintf(viewer,"\n");
1107: }
1108: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1109: PetscViewerASCIIPrintf(viewer,"];\n");
1110: }
1111: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1112: }
1113: PetscViewerFlush(viewer);
1114: return(0);
1115: }
1119: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1120: {
1121: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1122: PetscErrorCode ierr;
1123: int fd;
1124: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1125: PetscScalar *v,*anonz,*vals;
1126: PetscViewerFormat format;
1129: PetscViewerBinaryGetDescriptor(viewer,&fd);
1131: PetscViewerGetFormat(viewer,&format);
1132: if (format == PETSC_VIEWER_NATIVE) {
1133: /* store the matrix as a dense matrix */
1134: PetscMalloc1(4,&col_lens);
1136: col_lens[0] = MAT_FILE_CLASSID;
1137: col_lens[1] = m;
1138: col_lens[2] = n;
1139: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1141: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1142: PetscFree(col_lens);
1144: /* write out matrix, by rows */
1145: PetscMalloc1(m*n+1,&vals);
1146: v = a->v;
1147: for (j=0; j<n; j++) {
1148: for (i=0; i<m; i++) {
1149: vals[j + i*n] = *v++;
1150: }
1151: }
1152: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1153: PetscFree(vals);
1154: } else {
1155: PetscMalloc1(4+nz,&col_lens);
1157: col_lens[0] = MAT_FILE_CLASSID;
1158: col_lens[1] = m;
1159: col_lens[2] = n;
1160: col_lens[3] = nz;
1162: /* store lengths of each row and write (including header) to file */
1163: for (i=0; i<m; i++) col_lens[4+i] = n;
1164: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1166: /* Possibly should write in smaller increments, not whole matrix at once? */
1167: /* store column indices (zero start index) */
1168: ict = 0;
1169: for (i=0; i<m; i++) {
1170: for (j=0; j<n; j++) col_lens[ict++] = j;
1171: }
1172: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1173: PetscFree(col_lens);
1175: /* store nonzero values */
1176: PetscMalloc1(nz+1,&anonz);
1177: ict = 0;
1178: for (i=0; i<m; i++) {
1179: v = a->v + i;
1180: for (j=0; j<n; j++) {
1181: anonz[ict++] = *v; v += a->lda;
1182: }
1183: }
1184: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1185: PetscFree(anonz);
1186: }
1187: return(0);
1188: }
1190: #include <petscdraw.h>
1193: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1194: {
1195: Mat A = (Mat) Aa;
1196: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1197: PetscErrorCode ierr;
1198: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1199: int color = PETSC_DRAW_WHITE;
1200: PetscScalar *v = a->v;
1201: PetscViewer viewer;
1202: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1203: PetscViewerFormat format;
1206: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1207: PetscViewerGetFormat(viewer,&format);
1208: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1210: /* Loop over matrix elements drawing boxes */
1212: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1213: PetscDrawCollectiveBegin(draw);
1214: /* Blue for negative and Red for positive */
1215: for (j = 0; j < n; j++) {
1216: x_l = j; x_r = x_l + 1.0;
1217: for (i = 0; i < m; i++) {
1218: y_l = m - i - 1.0;
1219: y_r = y_l + 1.0;
1220: if (PetscRealPart(v[j*m+i]) > 0.) {
1221: color = PETSC_DRAW_RED;
1222: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1223: color = PETSC_DRAW_BLUE;
1224: } else {
1225: continue;
1226: }
1227: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1228: }
1229: }
1230: PetscDrawCollectiveEnd(draw);
1231: } else {
1232: /* use contour shading to indicate magnitude of values */
1233: /* first determine max of all nonzero values */
1234: PetscReal minv = 0.0, maxv = 0.0;
1235: PetscDraw popup;
1237: for (i=0; i < m*n; i++) {
1238: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1239: }
1240: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1241: PetscDrawGetPopup(draw,&popup);
1242: PetscDrawScalePopup(popup,minv,maxv);
1244: PetscDrawCollectiveBegin(draw);
1245: for (j=0; j<n; j++) {
1246: x_l = j;
1247: x_r = x_l + 1.0;
1248: for (i=0; i<m; i++) {
1249: y_l = m - i - 1.0;
1250: y_r = y_l + 1.0;
1251: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1252: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1253: }
1254: }
1255: PetscDrawCollectiveEnd(draw);
1256: }
1257: return(0);
1258: }
1262: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1263: {
1264: PetscDraw draw;
1265: PetscBool isnull;
1266: PetscReal xr,yr,xl,yl,h,w;
1270: PetscViewerDrawGetDraw(viewer,0,&draw);
1271: PetscDrawIsNull(draw,&isnull);
1272: if (isnull) return(0);
1274: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1275: xr += w; yr += h; xl = -w; yl = -h;
1276: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1277: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1278: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1279: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1280: PetscDrawSave(draw);
1281: return(0);
1282: }
1286: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1287: {
1289: PetscBool iascii,isbinary,isdraw;
1292: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1293: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1294: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1296: if (iascii) {
1297: MatView_SeqDense_ASCII(A,viewer);
1298: } else if (isbinary) {
1299: MatView_SeqDense_Binary(A,viewer);
1300: } else if (isdraw) {
1301: MatView_SeqDense_Draw(A,viewer);
1302: }
1303: return(0);
1304: }
1308: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1309: {
1310: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1314: #if defined(PETSC_USE_LOG)
1315: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1316: #endif
1317: PetscFree(l->pivots);
1318: MatDestroy(&l->ptapwork);
1319: if (!l->user_alloc) {PetscFree(l->v);}
1320: PetscFree(mat->data);
1322: PetscObjectChangeTypeName((PetscObject)mat,0);
1323: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1324: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1325: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1326: #if defined(PETSC_HAVE_ELEMENTAL)
1327: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1328: #endif
1329: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1330: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1331: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1332: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1333: PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1334: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1335: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1336: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1337: return(0);
1338: }
1342: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1343: {
1344: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1346: PetscInt k,j,m,n,M;
1347: PetscScalar *v,tmp;
1350: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1351: if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1352: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1353: else {
1354: for (j=0; j<m; j++) {
1355: for (k=0; k<j; k++) {
1356: tmp = v[j + k*M];
1357: v[j + k*M] = v[k + j*M];
1358: v[k + j*M] = tmp;
1359: }
1360: }
1361: }
1362: } else { /* out-of-place transpose */
1363: Mat tmat;
1364: Mat_SeqDense *tmatd;
1365: PetscScalar *v2;
1366: PetscInt M2;
1368: if (reuse == MAT_INITIAL_MATRIX) {
1369: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1370: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1371: MatSetType(tmat,((PetscObject)A)->type_name);
1372: MatSeqDenseSetPreallocation(tmat,NULL);
1373: } else {
1374: tmat = *matout;
1375: }
1376: tmatd = (Mat_SeqDense*)tmat->data;
1377: v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1378: for (j=0; j<n; j++) {
1379: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1380: }
1381: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1382: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1383: *matout = tmat;
1384: }
1385: return(0);
1386: }
1390: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1391: {
1392: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1393: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1394: PetscInt i,j;
1395: PetscScalar *v1,*v2;
1398: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1399: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1400: for (i=0; i<A1->rmap->n; i++) {
1401: v1 = mat1->v+i; v2 = mat2->v+i;
1402: for (j=0; j<A1->cmap->n; j++) {
1403: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1404: v1 += mat1->lda; v2 += mat2->lda;
1405: }
1406: }
1407: *flg = PETSC_TRUE;
1408: return(0);
1409: }
1413: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1414: {
1415: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1417: PetscInt i,n,len;
1418: PetscScalar *x,zero = 0.0;
1421: VecSet(v,zero);
1422: VecGetSize(v,&n);
1423: VecGetArray(v,&x);
1424: len = PetscMin(A->rmap->n,A->cmap->n);
1425: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1426: for (i=0; i<len; i++) {
1427: x[i] = mat->v[i*mat->lda + i];
1428: }
1429: VecRestoreArray(v,&x);
1430: return(0);
1431: }
1435: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1436: {
1437: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1438: const PetscScalar *l,*r;
1439: PetscScalar x,*v;
1440: PetscErrorCode ierr;
1441: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1444: if (ll) {
1445: VecGetSize(ll,&m);
1446: VecGetArrayRead(ll,&l);
1447: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1448: for (i=0; i<m; i++) {
1449: x = l[i];
1450: v = mat->v + i;
1451: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1452: }
1453: VecRestoreArrayRead(ll,&l);
1454: PetscLogFlops(1.0*n*m);
1455: }
1456: if (rr) {
1457: VecGetSize(rr,&n);
1458: VecGetArrayRead(rr,&r);
1459: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1460: for (i=0; i<n; i++) {
1461: x = r[i];
1462: v = mat->v + i*mat->lda;
1463: for (j=0; j<m; j++) (*v++) *= x;
1464: }
1465: VecRestoreArrayRead(rr,&r);
1466: PetscLogFlops(1.0*n*m);
1467: }
1468: return(0);
1469: }
1473: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1474: {
1475: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1476: PetscScalar *v = mat->v;
1477: PetscReal sum = 0.0;
1478: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1482: if (type == NORM_FROBENIUS) {
1483: if (lda>m) {
1484: for (j=0; j<A->cmap->n; j++) {
1485: v = mat->v+j*lda;
1486: for (i=0; i<m; i++) {
1487: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1488: }
1489: }
1490: } else {
1491: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1492: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1493: }
1494: }
1495: *nrm = PetscSqrtReal(sum);
1496: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1497: } else if (type == NORM_1) {
1498: *nrm = 0.0;
1499: for (j=0; j<A->cmap->n; j++) {
1500: v = mat->v + j*mat->lda;
1501: sum = 0.0;
1502: for (i=0; i<A->rmap->n; i++) {
1503: sum += PetscAbsScalar(*v); v++;
1504: }
1505: if (sum > *nrm) *nrm = sum;
1506: }
1507: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1508: } else if (type == NORM_INFINITY) {
1509: *nrm = 0.0;
1510: for (j=0; j<A->rmap->n; j++) {
1511: v = mat->v + j;
1512: sum = 0.0;
1513: for (i=0; i<A->cmap->n; i++) {
1514: sum += PetscAbsScalar(*v); v += mat->lda;
1515: }
1516: if (sum > *nrm) *nrm = sum;
1517: }
1518: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1519: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1520: return(0);
1521: }
1525: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1526: {
1527: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1531: switch (op) {
1532: case MAT_ROW_ORIENTED:
1533: aij->roworiented = flg;
1534: break;
1535: case MAT_NEW_NONZERO_LOCATIONS:
1536: case MAT_NEW_NONZERO_LOCATION_ERR:
1537: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1538: case MAT_NEW_DIAGONALS:
1539: case MAT_KEEP_NONZERO_PATTERN:
1540: case MAT_IGNORE_OFF_PROC_ENTRIES:
1541: case MAT_USE_HASH_TABLE:
1542: case MAT_IGNORE_LOWER_TRIANGULAR:
1543: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1544: break;
1545: case MAT_SPD:
1546: case MAT_SYMMETRIC:
1547: case MAT_STRUCTURALLY_SYMMETRIC:
1548: case MAT_HERMITIAN:
1549: case MAT_SYMMETRY_ETERNAL:
1550: /* These options are handled directly by MatSetOption() */
1551: break;
1552: default:
1553: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1554: }
1555: return(0);
1556: }
1560: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1561: {
1562: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1564: PetscInt lda=l->lda,m=A->rmap->n,j;
1567: if (lda>m) {
1568: for (j=0; j<A->cmap->n; j++) {
1569: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1570: }
1571: } else {
1572: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1573: }
1574: return(0);
1575: }
1579: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1580: {
1581: PetscErrorCode ierr;
1582: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1583: PetscInt m = l->lda, n = A->cmap->n, i,j;
1584: PetscScalar *slot,*bb;
1585: const PetscScalar *xx;
1588: #if defined(PETSC_USE_DEBUG)
1589: for (i=0; i<N; i++) {
1590: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1591: 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);
1592: }
1593: #endif
1595: /* fix right hand side if needed */
1596: if (x && b) {
1597: VecGetArrayRead(x,&xx);
1598: VecGetArray(b,&bb);
1599: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1600: VecRestoreArrayRead(x,&xx);
1601: VecRestoreArray(b,&bb);
1602: }
1604: for (i=0; i<N; i++) {
1605: slot = l->v + rows[i];
1606: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1607: }
1608: if (diag != 0.0) {
1609: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1610: for (i=0; i<N; i++) {
1611: slot = l->v + (m+1)*rows[i];
1612: *slot = diag;
1613: }
1614: }
1615: return(0);
1616: }
1620: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1621: {
1622: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1625: if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows");
1626: *array = mat->v;
1627: return(0);
1628: }
1632: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1633: {
1635: *array = 0; /* user cannot accidently use the array later */
1636: return(0);
1637: }
1641: /*@C
1642: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1644: Not Collective
1646: Input Parameter:
1647: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1649: Output Parameter:
1650: . array - pointer to the data
1652: Level: intermediate
1654: .seealso: MatDenseRestoreArray()
1655: @*/
1656: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1657: {
1661: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1662: return(0);
1663: }
1667: /*@C
1668: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1670: Not Collective
1672: Input Parameters:
1673: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1674: . array - pointer to the data
1676: Level: intermediate
1678: .seealso: MatDenseGetArray()
1679: @*/
1680: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1681: {
1685: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1686: return(0);
1687: }
1691: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1692: {
1693: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1695: PetscInt i,j,nrows,ncols;
1696: const PetscInt *irow,*icol;
1697: PetscScalar *av,*bv,*v = mat->v;
1698: Mat newmat;
1701: ISGetIndices(isrow,&irow);
1702: ISGetIndices(iscol,&icol);
1703: ISGetLocalSize(isrow,&nrows);
1704: ISGetLocalSize(iscol,&ncols);
1706: /* Check submatrixcall */
1707: if (scall == MAT_REUSE_MATRIX) {
1708: PetscInt n_cols,n_rows;
1709: MatGetSize(*B,&n_rows,&n_cols);
1710: if (n_rows != nrows || n_cols != ncols) {
1711: /* resize the result matrix to match number of requested rows/columns */
1712: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1713: }
1714: newmat = *B;
1715: } else {
1716: /* Create and fill new matrix */
1717: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1718: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1719: MatSetType(newmat,((PetscObject)A)->type_name);
1720: MatSeqDenseSetPreallocation(newmat,NULL);
1721: }
1723: /* Now extract the data pointers and do the copy,column at a time */
1724: bv = ((Mat_SeqDense*)newmat->data)->v;
1726: for (i=0; i<ncols; i++) {
1727: av = v + mat->lda*icol[i];
1728: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1729: }
1731: /* Assemble the matrices so that the correct flags are set */
1732: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1733: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1735: /* Free work space */
1736: ISRestoreIndices(isrow,&irow);
1737: ISRestoreIndices(iscol,&icol);
1738: *B = newmat;
1739: return(0);
1740: }
1744: static PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1745: {
1747: PetscInt i;
1750: if (scall == MAT_INITIAL_MATRIX) {
1751: PetscMalloc1(n+1,B);
1752: }
1754: for (i=0; i<n; i++) {
1755: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1756: }
1757: return(0);
1758: }
1762: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1763: {
1765: return(0);
1766: }
1770: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1771: {
1773: return(0);
1774: }
1778: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1779: {
1780: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1782: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1785: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1786: if (A->ops->copy != B->ops->copy) {
1787: MatCopy_Basic(A,B,str);
1788: return(0);
1789: }
1790: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1791: if (lda1>m || lda2>m) {
1792: for (j=0; j<n; j++) {
1793: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1794: }
1795: } else {
1796: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1797: }
1798: return(0);
1799: }
1803: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1804: {
1808: MatSeqDenseSetPreallocation(A,0);
1809: return(0);
1810: }
1814: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1815: {
1816: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1817: PetscInt i,nz = A->rmap->n*A->cmap->n;
1818: PetscScalar *aa = a->v;
1821: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1822: return(0);
1823: }
1827: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1828: {
1829: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1830: PetscInt i,nz = A->rmap->n*A->cmap->n;
1831: PetscScalar *aa = a->v;
1834: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1835: return(0);
1836: }
1840: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1841: {
1842: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1843: PetscInt i,nz = A->rmap->n*A->cmap->n;
1844: PetscScalar *aa = a->v;
1847: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1848: return(0);
1849: }
1851: /* ----------------------------------------------------------------*/
1854: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1855: {
1859: if (scall == MAT_INITIAL_MATRIX) {
1860: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1861: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1862: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1863: }
1864: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1865: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1866: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1867: return(0);
1868: }
1872: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1873: {
1875: PetscInt m=A->rmap->n,n=B->cmap->n;
1876: Mat Cmat;
1879: if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n);
1880: MatCreate(PETSC_COMM_SELF,&Cmat);
1881: MatSetSizes(Cmat,m,n,m,n);
1882: MatSetType(Cmat,MATSEQDENSE);
1883: MatSeqDenseSetPreallocation(Cmat,NULL);
1885: *C = Cmat;
1886: return(0);
1887: }
1891: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1892: {
1893: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1894: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1895: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1896: PetscBLASInt m,n,k;
1897: PetscScalar _DOne=1.0,_DZero=0.0;
1899: PetscBool flg;
1902: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
1903: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense");
1905: /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1906: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1907: if (flg) {
1908: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1909: (*C->ops->matmultnumeric)(A,B,C);
1910: return(0);
1911: }
1913: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1914: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1915: PetscBLASIntCast(A->rmap->n,&m);
1916: PetscBLASIntCast(B->cmap->n,&n);
1917: PetscBLASIntCast(A->cmap->n,&k);
1918: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1919: return(0);
1920: }
1924: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1925: {
1929: if (scall == MAT_INITIAL_MATRIX) {
1930: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1931: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1932: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1933: }
1934: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1935: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1936: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1937: return(0);
1938: }
1942: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1943: {
1945: PetscInt m=A->cmap->n,n=B->cmap->n;
1946: Mat Cmat;
1949: if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n);
1950: MatCreate(PETSC_COMM_SELF,&Cmat);
1951: MatSetSizes(Cmat,m,n,m,n);
1952: MatSetType(Cmat,MATSEQDENSE);
1953: MatSeqDenseSetPreallocation(Cmat,NULL);
1955: Cmat->assembled = PETSC_TRUE;
1957: *C = Cmat;
1958: return(0);
1959: }
1963: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1964: {
1965: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1966: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1967: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1968: PetscBLASInt m,n,k;
1969: PetscScalar _DOne=1.0,_DZero=0.0;
1973: PetscBLASIntCast(A->cmap->n,&m);
1974: PetscBLASIntCast(B->cmap->n,&n);
1975: PetscBLASIntCast(A->rmap->n,&k);
1976: /*
1977: Note the m and n arguments below are the number rows and columns of A', not A!
1978: */
1979: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1980: return(0);
1981: }
1985: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1986: {
1987: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1989: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1990: PetscScalar *x;
1991: MatScalar *aa = a->v;
1994: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1996: VecSet(v,0.0);
1997: VecGetArray(v,&x);
1998: VecGetLocalSize(v,&p);
1999: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2000: for (i=0; i<m; i++) {
2001: x[i] = aa[i]; if (idx) idx[i] = 0;
2002: for (j=1; j<n; j++) {
2003: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2004: }
2005: }
2006: VecRestoreArray(v,&x);
2007: return(0);
2008: }
2012: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2013: {
2014: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2016: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2017: PetscScalar *x;
2018: PetscReal atmp;
2019: MatScalar *aa = a->v;
2022: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2024: VecSet(v,0.0);
2025: VecGetArray(v,&x);
2026: VecGetLocalSize(v,&p);
2027: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2028: for (i=0; i<m; i++) {
2029: x[i] = PetscAbsScalar(aa[i]);
2030: for (j=1; j<n; j++) {
2031: atmp = PetscAbsScalar(aa[i+m*j]);
2032: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2033: }
2034: }
2035: VecRestoreArray(v,&x);
2036: return(0);
2037: }
2041: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2042: {
2043: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2045: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2046: PetscScalar *x;
2047: MatScalar *aa = a->v;
2050: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2052: VecSet(v,0.0);
2053: VecGetArray(v,&x);
2054: VecGetLocalSize(v,&p);
2055: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2056: for (i=0; i<m; i++) {
2057: x[i] = aa[i]; if (idx) idx[i] = 0;
2058: for (j=1; j<n; j++) {
2059: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2060: }
2061: }
2062: VecRestoreArray(v,&x);
2063: return(0);
2064: }
2068: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2069: {
2070: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2072: PetscScalar *x;
2075: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2077: VecGetArray(v,&x);
2078: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2079: VecRestoreArray(v,&x);
2080: return(0);
2081: }
2086: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2087: {
2089: PetscInt i,j,m,n;
2090: PetscScalar *a;
2093: MatGetSize(A,&m,&n);
2094: PetscMemzero(norms,n*sizeof(PetscReal));
2095: MatDenseGetArray(A,&a);
2096: if (type == NORM_2) {
2097: for (i=0; i<n; i++) {
2098: for (j=0; j<m; j++) {
2099: norms[i] += PetscAbsScalar(a[j]*a[j]);
2100: }
2101: a += m;
2102: }
2103: } else if (type == NORM_1) {
2104: for (i=0; i<n; i++) {
2105: for (j=0; j<m; j++) {
2106: norms[i] += PetscAbsScalar(a[j]);
2107: }
2108: a += m;
2109: }
2110: } else if (type == NORM_INFINITY) {
2111: for (i=0; i<n; i++) {
2112: for (j=0; j<m; j++) {
2113: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2114: }
2115: a += m;
2116: }
2117: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2118: MatDenseRestoreArray(A,&a);
2119: if (type == NORM_2) {
2120: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2121: }
2122: return(0);
2123: }
2127: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2128: {
2130: PetscScalar *a;
2131: PetscInt m,n,i;
2134: MatGetSize(x,&m,&n);
2135: MatDenseGetArray(x,&a);
2136: for (i=0; i<m*n; i++) {
2137: PetscRandomGetValue(rctx,a+i);
2138: }
2139: MatDenseRestoreArray(x,&a);
2140: return(0);
2141: }
2145: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2146: {
2148: *missing = PETSC_FALSE;
2149: return(0);
2150: }
2153: /* -------------------------------------------------------------------*/
2154: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2155: MatGetRow_SeqDense,
2156: MatRestoreRow_SeqDense,
2157: MatMult_SeqDense,
2158: /* 4*/ MatMultAdd_SeqDense,
2159: MatMultTranspose_SeqDense,
2160: MatMultTransposeAdd_SeqDense,
2161: 0,
2162: 0,
2163: 0,
2164: /* 10*/ 0,
2165: MatLUFactor_SeqDense,
2166: MatCholeskyFactor_SeqDense,
2167: MatSOR_SeqDense,
2168: MatTranspose_SeqDense,
2169: /* 15*/ MatGetInfo_SeqDense,
2170: MatEqual_SeqDense,
2171: MatGetDiagonal_SeqDense,
2172: MatDiagonalScale_SeqDense,
2173: MatNorm_SeqDense,
2174: /* 20*/ MatAssemblyBegin_SeqDense,
2175: MatAssemblyEnd_SeqDense,
2176: MatSetOption_SeqDense,
2177: MatZeroEntries_SeqDense,
2178: /* 24*/ MatZeroRows_SeqDense,
2179: 0,
2180: 0,
2181: 0,
2182: 0,
2183: /* 29*/ MatSetUp_SeqDense,
2184: 0,
2185: 0,
2186: 0,
2187: 0,
2188: /* 34*/ MatDuplicate_SeqDense,
2189: 0,
2190: 0,
2191: 0,
2192: 0,
2193: /* 39*/ MatAXPY_SeqDense,
2194: MatGetSubMatrices_SeqDense,
2195: 0,
2196: MatGetValues_SeqDense,
2197: MatCopy_SeqDense,
2198: /* 44*/ MatGetRowMax_SeqDense,
2199: MatScale_SeqDense,
2200: MatShift_Basic,
2201: 0,
2202: MatZeroRowsColumns_SeqDense,
2203: /* 49*/ MatSetRandom_SeqDense,
2204: 0,
2205: 0,
2206: 0,
2207: 0,
2208: /* 54*/ 0,
2209: 0,
2210: 0,
2211: 0,
2212: 0,
2213: /* 59*/ 0,
2214: MatDestroy_SeqDense,
2215: MatView_SeqDense,
2216: 0,
2217: 0,
2218: /* 64*/ 0,
2219: 0,
2220: 0,
2221: 0,
2222: 0,
2223: /* 69*/ MatGetRowMaxAbs_SeqDense,
2224: 0,
2225: 0,
2226: 0,
2227: 0,
2228: /* 74*/ 0,
2229: 0,
2230: 0,
2231: 0,
2232: 0,
2233: /* 79*/ 0,
2234: 0,
2235: 0,
2236: 0,
2237: /* 83*/ MatLoad_SeqDense,
2238: 0,
2239: MatIsHermitian_SeqDense,
2240: 0,
2241: 0,
2242: 0,
2243: /* 89*/ MatMatMult_SeqDense_SeqDense,
2244: MatMatMultSymbolic_SeqDense_SeqDense,
2245: MatMatMultNumeric_SeqDense_SeqDense,
2246: MatPtAP_SeqDense_SeqDense,
2247: MatPtAPSymbolic_SeqDense_SeqDense,
2248: /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2249: 0,
2250: 0,
2251: 0,
2252: 0,
2253: /* 99*/ 0,
2254: 0,
2255: 0,
2256: MatConjugate_SeqDense,
2257: 0,
2258: /*104*/ 0,
2259: MatRealPart_SeqDense,
2260: MatImaginaryPart_SeqDense,
2261: 0,
2262: 0,
2263: /*109*/ MatMatSolve_SeqDense,
2264: 0,
2265: MatGetRowMin_SeqDense,
2266: MatGetColumnVector_SeqDense,
2267: MatMissingDiagonal_SeqDense,
2268: /*114*/ 0,
2269: 0,
2270: 0,
2271: 0,
2272: 0,
2273: /*119*/ 0,
2274: 0,
2275: 0,
2276: 0,
2277: 0,
2278: /*124*/ 0,
2279: MatGetColumnNorms_SeqDense,
2280: 0,
2281: 0,
2282: 0,
2283: /*129*/ 0,
2284: MatTransposeMatMult_SeqDense_SeqDense,
2285: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2286: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2287: 0,
2288: /*134*/ 0,
2289: 0,
2290: 0,
2291: 0,
2292: 0,
2293: /*139*/ 0,
2294: 0,
2295: 0
2296: };
2300: /*@C
2301: MatCreateSeqDense - Creates a sequential dense matrix that
2302: is stored in column major order (the usual Fortran 77 manner). Many
2303: of the matrix operations use the BLAS and LAPACK routines.
2305: Collective on MPI_Comm
2307: Input Parameters:
2308: + comm - MPI communicator, set to PETSC_COMM_SELF
2309: . m - number of rows
2310: . n - number of columns
2311: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2312: to control all matrix memory allocation.
2314: Output Parameter:
2315: . A - the matrix
2317: Notes:
2318: The data input variable is intended primarily for Fortran programmers
2319: who wish to allocate their own matrix memory space. Most users should
2320: set data=NULL.
2322: Level: intermediate
2324: .keywords: dense, matrix, LAPACK, BLAS
2326: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2327: @*/
2328: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2329: {
2333: MatCreate(comm,A);
2334: MatSetSizes(*A,m,n,m,n);
2335: MatSetType(*A,MATSEQDENSE);
2336: MatSeqDenseSetPreallocation(*A,data);
2337: return(0);
2338: }
2342: /*@C
2343: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2345: Collective on MPI_Comm
2347: Input Parameters:
2348: + B - the matrix
2349: - data - the array (or NULL)
2351: Notes:
2352: The data input variable is intended primarily for Fortran programmers
2353: who wish to allocate their own matrix memory space. Most users should
2354: need not call this routine.
2356: Level: intermediate
2358: .keywords: dense, matrix, LAPACK, BLAS
2360: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2362: @*/
2363: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2364: {
2368: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2369: return(0);
2370: }
2374: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2375: {
2376: Mat_SeqDense *b;
2380: B->preallocated = PETSC_TRUE;
2382: PetscLayoutSetUp(B->rmap);
2383: PetscLayoutSetUp(B->cmap);
2385: b = (Mat_SeqDense*)B->data;
2386: b->Mmax = B->rmap->n;
2387: b->Nmax = B->cmap->n;
2388: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2390: if (!data) { /* petsc-allocated storage */
2391: if (!b->user_alloc) { PetscFree(b->v); }
2392: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2393: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2395: b->user_alloc = PETSC_FALSE;
2396: } else { /* user-allocated storage */
2397: if (!b->user_alloc) { PetscFree(b->v); }
2398: b->v = data;
2399: b->user_alloc = PETSC_TRUE;
2400: }
2401: B->assembled = PETSC_TRUE;
2402: return(0);
2403: }
2405: #if defined(PETSC_HAVE_ELEMENTAL)
2408: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2409: {
2410: Mat mat_elemental;
2412: PetscScalar *array,*v_colwise;
2413: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2416: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2417: MatDenseGetArray(A,&array);
2418: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2419: k = 0;
2420: for (j=0; j<N; j++) {
2421: cols[j] = j;
2422: for (i=0; i<M; i++) {
2423: v_colwise[j*M+i] = array[k++];
2424: }
2425: }
2426: for (i=0; i<M; i++) {
2427: rows[i] = i;
2428: }
2429: MatDenseRestoreArray(A,&array);
2431: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2432: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2433: MatSetType(mat_elemental,MATELEMENTAL);
2434: MatSetUp(mat_elemental);
2436: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2437: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2438: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2439: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2440: PetscFree3(v_colwise,rows,cols);
2442: if (reuse == MAT_INPLACE_MATRIX) {
2443: MatHeaderReplace(A,&mat_elemental);
2444: } else {
2445: *newmat = mat_elemental;
2446: }
2447: return(0);
2448: }
2449: #endif
2453: /*@C
2454: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2456: Input parameter:
2457: + A - the matrix
2458: - lda - the leading dimension
2460: Notes:
2461: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2462: it asserts that the preallocation has a leading dimension (the LDA parameter
2463: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2465: Level: intermediate
2467: .keywords: dense, matrix, LAPACK, BLAS
2469: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2471: @*/
2472: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2473: {
2474: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2477: 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);
2478: b->lda = lda;
2479: b->changelda = PETSC_FALSE;
2480: b->Mmax = PetscMax(b->Mmax,lda);
2481: return(0);
2482: }
2484: /*MC
2485: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2487: Options Database Keys:
2488: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2490: Level: beginner
2492: .seealso: MatCreateSeqDense()
2494: M*/
2498: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2499: {
2500: Mat_SeqDense *b;
2502: PetscMPIInt size;
2505: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2506: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2508: PetscNewLog(B,&b);
2509: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2510: B->data = (void*)b;
2512: b->pivots = 0;
2513: b->roworiented = PETSC_TRUE;
2514: b->v = 0;
2515: b->changelda = PETSC_FALSE;
2517: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2518: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2519: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2520: #if defined(PETSC_HAVE_ELEMENTAL)
2521: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2522: #endif
2523: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2524: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2525: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2526: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2527: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2528: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2529: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2530: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2531: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2532: return(0);
2533: }