Actual source code: dense.c
petsc-3.6.1 2015-08-06
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: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
14: {
15: Mat B;
16: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
17: Mat_SeqDense *b;
19: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
20: MatScalar *av=a->a;
23: MatCreate(PetscObjectComm((PetscObject)A),&B);
24: MatSetSizes(B,m,n,m,n);
25: MatSetType(B,MATSEQDENSE);
26: MatSeqDenseSetPreallocation(B,NULL);
28: b = (Mat_SeqDense*)(B->data);
29: for (i=0; i<m; i++) {
30: PetscInt j;
31: for (j=0;j<ai[1]-ai[0];j++) {
32: b->v[*aj*m+i] = *av;
33: aj++;
34: av++;
35: }
36: ai++;
37: }
38: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
39: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
41: if (reuse == MAT_REUSE_MATRIX) {
42: MatHeaderReplace(A,B);
43: } else {
44: *newmat = B;
45: }
46: return(0);
47: }
51: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
52: {
53: Mat B;
54: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
56: PetscInt i, j;
57: PetscInt *rows, *nnz;
58: MatScalar *aa = a->v, *vals;
61: MatCreate(PetscObjectComm((PetscObject)A),&B);
62: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
63: MatSetType(B,MATSEQAIJ);
64: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
65: for (j=0; j<A->cmap->n; j++) {
66: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
67: aa += a->lda;
68: }
69: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
70: aa = a->v;
71: for (j=0; j<A->cmap->n; j++) {
72: PetscInt numRows = 0;
73: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
74: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
75: aa += a->lda;
76: }
77: PetscFree3(rows,nnz,vals);
78: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
81: if (reuse == MAT_REUSE_MATRIX) {
82: MatHeaderReplace(A,B);
83: } else {
84: *newmat = B;
85: }
86: return(0);
87: }
91: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
92: {
93: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
94: PetscScalar oalpha = alpha;
95: PetscInt j;
96: PetscBLASInt N,m,ldax,lday,one = 1;
100: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
101: PetscBLASIntCast(X->rmap->n,&m);
102: PetscBLASIntCast(x->lda,&ldax);
103: PetscBLASIntCast(y->lda,&lday);
104: if (ldax>m || lday>m) {
105: for (j=0; j<X->cmap->n; j++) {
106: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
107: }
108: } else {
109: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
110: }
111: PetscObjectStateIncrease((PetscObject)Y);
112: PetscLogFlops(PetscMax(2*N-1,0));
113: return(0);
114: }
118: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
119: {
120: PetscInt N = A->rmap->n*A->cmap->n;
123: info->block_size = 1.0;
124: info->nz_allocated = (double)N;
125: info->nz_used = (double)N;
126: info->nz_unneeded = (double)0;
127: info->assemblies = (double)A->num_ass;
128: info->mallocs = 0;
129: info->memory = ((PetscObject)A)->mem;
130: info->fill_ratio_given = 0;
131: info->fill_ratio_needed = 0;
132: info->factor_mallocs = 0;
133: return(0);
134: }
138: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
139: {
140: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
141: PetscScalar oalpha = alpha;
143: PetscBLASInt one = 1,j,nz,lda;
146: PetscBLASIntCast(a->lda,&lda);
147: if (lda>A->rmap->n) {
148: PetscBLASIntCast(A->rmap->n,&nz);
149: for (j=0; j<A->cmap->n; j++) {
150: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
151: }
152: } else {
153: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
154: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
155: }
156: PetscLogFlops(nz);
157: return(0);
158: }
162: PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
163: {
164: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
165: PetscInt i,j,m = A->rmap->n,N;
166: PetscScalar *v = a->v;
169: *fl = PETSC_FALSE;
170: if (A->rmap->n != A->cmap->n) return(0);
171: N = a->lda;
173: for (i=0; i<m; i++) {
174: for (j=i+1; j<m; j++) {
175: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
176: }
177: }
178: *fl = PETSC_TRUE;
179: return(0);
180: }
184: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
185: {
186: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
188: PetscInt lda = (PetscInt)mat->lda,j,m;
191: PetscLayoutReference(A->rmap,&newi->rmap);
192: PetscLayoutReference(A->cmap,&newi->cmap);
193: MatSeqDenseSetPreallocation(newi,NULL);
194: if (cpvalues == MAT_COPY_VALUES) {
195: l = (Mat_SeqDense*)newi->data;
196: if (lda>A->rmap->n) {
197: m = A->rmap->n;
198: for (j=0; j<A->cmap->n; j++) {
199: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
200: }
201: } else {
202: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
203: }
204: }
205: newi->assembled = PETSC_TRUE;
206: return(0);
207: }
211: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
212: {
216: MatCreate(PetscObjectComm((PetscObject)A),newmat);
217: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
218: MatSetType(*newmat,((PetscObject)A)->type_name);
219: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
220: return(0);
221: }
224: extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
228: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
229: {
230: MatFactorInfo info;
234: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
235: MatLUFactor_SeqDense(fact,0,0,&info);
236: return(0);
237: }
241: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
242: {
243: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
244: PetscErrorCode ierr;
245: const PetscScalar *x;
246: PetscScalar *y;
247: PetscBLASInt one = 1,info,m;
250: PetscBLASIntCast(A->rmap->n,&m);
251: VecGetArrayRead(xx,&x);
252: VecGetArray(yy,&y);
253: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
254: if (A->factortype == MAT_FACTOR_LU) {
255: #if defined(PETSC_MISSING_LAPACK_GETRS)
256: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
257: #else
258: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
259: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
260: #endif
261: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
262: #if defined(PETSC_MISSING_LAPACK_POTRS)
263: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
264: #else
265: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
266: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
267: #endif
268: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
269: VecRestoreArrayRead(xx,&x);
270: VecRestoreArray(yy,&y);
271: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
272: return(0);
273: }
277: PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
278: {
279: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
281: PetscScalar *b,*x;
282: PetscInt n;
283: PetscBLASInt nrhs,info,m;
284: PetscBool flg;
287: PetscBLASIntCast(A->rmap->n,&m);
288: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
289: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
290: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
291: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
293: MatGetSize(B,NULL,&n);
294: PetscBLASIntCast(n,&nrhs);
295: MatDenseGetArray(B,&b);
296: MatDenseGetArray(X,&x);
298: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
300: if (A->factortype == MAT_FACTOR_LU) {
301: #if defined(PETSC_MISSING_LAPACK_GETRS)
302: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
303: #else
304: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
305: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
306: #endif
307: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
308: #if defined(PETSC_MISSING_LAPACK_POTRS)
309: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
310: #else
311: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
312: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
313: #endif
314: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
316: MatDenseRestoreArray(B,&b);
317: MatDenseRestoreArray(X,&x);
318: PetscLogFlops(nrhs*(2.0*m*m - m));
319: return(0);
320: }
324: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
325: {
326: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
327: PetscErrorCode ierr;
328: const PetscScalar *x;
329: PetscScalar *y;
330: PetscBLASInt one = 1,info,m;
333: PetscBLASIntCast(A->rmap->n,&m);
334: VecGetArrayRead(xx,&x);
335: VecGetArray(yy,&y);
336: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
337: /* assume if pivots exist then use LU; else Cholesky */
338: if (mat->pivots) {
339: #if defined(PETSC_MISSING_LAPACK_GETRS)
340: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
341: #else
342: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
343: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
344: #endif
345: } else {
346: #if defined(PETSC_MISSING_LAPACK_POTRS)
347: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
348: #else
349: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
350: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
351: #endif
352: }
353: VecRestoreArrayRead(xx,&x);
354: VecRestoreArray(yy,&y);
355: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
356: return(0);
357: }
361: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
362: {
363: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
364: PetscErrorCode ierr;
365: const PetscScalar *x;
366: PetscScalar *y,sone = 1.0;
367: Vec tmp = 0;
368: PetscBLASInt one = 1,info,m;
371: PetscBLASIntCast(A->rmap->n,&m);
372: VecGetArrayRead(xx,&x);
373: VecGetArray(yy,&y);
374: if (!A->rmap->n || !A->cmap->n) return(0);
375: if (yy == zz) {
376: VecDuplicate(yy,&tmp);
377: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
378: VecCopy(yy,tmp);
379: }
380: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
381: /* assume if pivots exist then use LU; else Cholesky */
382: if (mat->pivots) {
383: #if defined(PETSC_MISSING_LAPACK_GETRS)
384: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
385: #else
386: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
387: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
388: #endif
389: } else {
390: #if defined(PETSC_MISSING_LAPACK_POTRS)
391: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
392: #else
393: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
394: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
395: #endif
396: }
397: if (tmp) {
398: VecAXPY(yy,sone,tmp);
399: VecDestroy(&tmp);
400: } else {
401: VecAXPY(yy,sone,zz);
402: }
403: VecRestoreArrayRead(xx,&x);
404: VecRestoreArray(yy,&y);
405: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
406: return(0);
407: }
411: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
412: {
413: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
414: PetscErrorCode ierr;
415: const PetscScalar *x;
416: PetscScalar *y,sone = 1.0;
417: Vec tmp;
418: PetscBLASInt one = 1,info,m;
421: PetscBLASIntCast(A->rmap->n,&m);
422: if (!A->rmap->n || !A->cmap->n) return(0);
423: VecGetArrayRead(xx,&x);
424: VecGetArray(yy,&y);
425: if (yy == zz) {
426: VecDuplicate(yy,&tmp);
427: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
428: VecCopy(yy,tmp);
429: }
430: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
431: /* assume if pivots exist then use LU; else Cholesky */
432: if (mat->pivots) {
433: #if defined(PETSC_MISSING_LAPACK_GETRS)
434: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
435: #else
436: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
437: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
438: #endif
439: } else {
440: #if defined(PETSC_MISSING_LAPACK_POTRS)
441: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
442: #else
443: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
444: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
445: #endif
446: }
447: if (tmp) {
448: VecAXPY(yy,sone,tmp);
449: VecDestroy(&tmp);
450: } else {
451: VecAXPY(yy,sone,zz);
452: }
453: VecRestoreArrayRead(xx,&x);
454: VecRestoreArray(yy,&y);
455: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
456: return(0);
457: }
459: /* ---------------------------------------------------------------*/
460: /* COMMENT: I have chosen to hide row permutation in the pivots,
461: rather than put it in the Mat->row slot.*/
464: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
465: {
466: #if defined(PETSC_MISSING_LAPACK_GETRF)
468: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
469: #else
470: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
472: PetscBLASInt n,m,info;
475: PetscBLASIntCast(A->cmap->n,&n);
476: PetscBLASIntCast(A->rmap->n,&m);
477: if (!mat->pivots) {
478: PetscMalloc1(A->rmap->n+1,&mat->pivots);
479: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
480: }
481: if (!A->rmap->n || !A->cmap->n) return(0);
482: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
484: PetscFPTrapPop();
486: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
487: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
488: A->ops->solve = MatSolve_SeqDense;
489: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
490: A->ops->solveadd = MatSolveAdd_SeqDense;
491: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
492: A->factortype = MAT_FACTOR_LU;
494: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
495: #endif
496: return(0);
497: }
501: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
502: {
503: #if defined(PETSC_MISSING_LAPACK_POTRF)
505: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
506: #else
507: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
509: PetscBLASInt info,n;
512: PetscBLASIntCast(A->cmap->n,&n);
513: PetscFree(mat->pivots);
515: if (!A->rmap->n || !A->cmap->n) return(0);
516: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
517: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
518: A->ops->solve = MatSolve_SeqDense;
519: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
520: A->ops->solveadd = MatSolveAdd_SeqDense;
521: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
522: A->factortype = MAT_FACTOR_CHOLESKY;
524: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
525: #endif
526: return(0);
527: }
532: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
533: {
535: MatFactorInfo info;
538: info.fill = 1.0;
540: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
541: MatCholeskyFactor_SeqDense(fact,0,&info);
542: return(0);
543: }
547: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
548: {
550: fact->assembled = PETSC_TRUE;
551: fact->preallocated = PETSC_TRUE;
552: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
553: return(0);
554: }
558: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
559: {
561: fact->preallocated = PETSC_TRUE;
562: fact->assembled = PETSC_TRUE;
563: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
564: return(0);
565: }
569: PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
570: {
574: MatCreate(PetscObjectComm((PetscObject)A),fact);
575: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
576: MatSetType(*fact,((PetscObject)A)->type_name);
577: if (ftype == MAT_FACTOR_LU) {
578: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
579: } else {
580: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
581: }
582: (*fact)->factortype = ftype;
583: return(0);
584: }
586: /* ------------------------------------------------------------------*/
589: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
590: {
591: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
592: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
593: const PetscScalar *b;
594: PetscErrorCode ierr;
595: PetscInt m = A->rmap->n,i;
596: PetscBLASInt o = 1,bm;
599: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
600: PetscBLASIntCast(m,&bm);
601: if (flag & SOR_ZERO_INITIAL_GUESS) {
602: /* this is a hack fix, should have another version without the second BLASdot */
603: VecSet(xx,zero);
604: }
605: VecGetArray(xx,&x);
606: VecGetArrayRead(bb,&b);
607: its = its*lits;
608: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
609: while (its--) {
610: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
611: for (i=0; i<m; i++) {
612: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
613: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
614: }
615: }
616: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
617: for (i=m-1; i>=0; i--) {
618: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
619: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
620: }
621: }
622: }
623: VecRestoreArrayRead(bb,&b);
624: VecRestoreArray(xx,&x);
625: return(0);
626: }
628: /* -----------------------------------------------------------------*/
631: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
632: {
633: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
634: const PetscScalar *v = mat->v,*x;
635: PetscScalar *y;
636: PetscErrorCode ierr;
637: PetscBLASInt m, n,_One=1;
638: PetscScalar _DOne=1.0,_DZero=0.0;
641: PetscBLASIntCast(A->rmap->n,&m);
642: PetscBLASIntCast(A->cmap->n,&n);
643: if (!A->rmap->n || !A->cmap->n) return(0);
644: VecGetArrayRead(xx,&x);
645: VecGetArray(yy,&y);
646: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
647: VecRestoreArrayRead(xx,&x);
648: VecRestoreArray(yy,&y);
649: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
650: return(0);
651: }
655: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
656: {
657: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
658: PetscScalar *y,_DOne=1.0,_DZero=0.0;
659: PetscErrorCode ierr;
660: PetscBLASInt m, n, _One=1;
661: const PetscScalar *v = mat->v,*x;
664: PetscBLASIntCast(A->rmap->n,&m);
665: PetscBLASIntCast(A->cmap->n,&n);
666: if (!A->rmap->n || !A->cmap->n) return(0);
667: VecGetArrayRead(xx,&x);
668: VecGetArray(yy,&y);
669: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
670: VecRestoreArrayRead(xx,&x);
671: VecRestoreArray(yy,&y);
672: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
673: return(0);
674: }
678: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
679: {
680: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
681: const PetscScalar *v = mat->v,*x;
682: PetscScalar *y,_DOne=1.0;
683: PetscErrorCode ierr;
684: PetscBLASInt m, n, _One=1;
687: PetscBLASIntCast(A->rmap->n,&m);
688: PetscBLASIntCast(A->cmap->n,&n);
689: if (!A->rmap->n || !A->cmap->n) return(0);
690: if (zz != yy) {VecCopy(zz,yy);}
691: VecGetArrayRead(xx,&x);
692: VecGetArray(yy,&y);
693: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
694: VecRestoreArrayRead(xx,&x);
695: VecRestoreArray(yy,&y);
696: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
697: return(0);
698: }
702: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
703: {
704: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
705: const PetscScalar *v = mat->v,*x;
706: PetscScalar *y;
707: PetscErrorCode ierr;
708: PetscBLASInt m, n, _One=1;
709: PetscScalar _DOne=1.0;
712: PetscBLASIntCast(A->rmap->n,&m);
713: PetscBLASIntCast(A->cmap->n,&n);
714: if (!A->rmap->n || !A->cmap->n) return(0);
715: if (zz != yy) {VecCopy(zz,yy);}
716: VecGetArrayRead(xx,&x);
717: VecGetArray(yy,&y);
718: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
719: VecRestoreArrayRead(xx,&x);
720: VecRestoreArray(yy,&y);
721: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
722: return(0);
723: }
725: /* -----------------------------------------------------------------*/
728: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
729: {
730: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
731: PetscScalar *v;
733: PetscInt i;
736: *ncols = A->cmap->n;
737: if (cols) {
738: PetscMalloc1(A->cmap->n+1,cols);
739: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
740: }
741: if (vals) {
742: PetscMalloc1(A->cmap->n+1,vals);
743: v = mat->v + row;
744: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
745: }
746: return(0);
747: }
751: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
752: {
756: if (cols) {PetscFree(*cols);}
757: if (vals) {PetscFree(*vals); }
758: return(0);
759: }
760: /* ----------------------------------------------------------------*/
763: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
764: {
765: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
766: PetscInt i,j,idx=0;
769: if (!mat->roworiented) {
770: if (addv == INSERT_VALUES) {
771: for (j=0; j<n; j++) {
772: if (indexn[j] < 0) {idx += m; continue;}
773: #if defined(PETSC_USE_DEBUG)
774: 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);
775: #endif
776: for (i=0; i<m; i++) {
777: if (indexm[i] < 0) {idx++; continue;}
778: #if defined(PETSC_USE_DEBUG)
779: 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);
780: #endif
781: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
782: }
783: }
784: } else {
785: for (j=0; j<n; j++) {
786: if (indexn[j] < 0) {idx += m; continue;}
787: #if defined(PETSC_USE_DEBUG)
788: 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);
789: #endif
790: for (i=0; i<m; i++) {
791: if (indexm[i] < 0) {idx++; continue;}
792: #if defined(PETSC_USE_DEBUG)
793: 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);
794: #endif
795: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
796: }
797: }
798: }
799: } else {
800: if (addv == INSERT_VALUES) {
801: for (i=0; i<m; i++) {
802: if (indexm[i] < 0) { idx += n; continue;}
803: #if defined(PETSC_USE_DEBUG)
804: 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);
805: #endif
806: for (j=0; j<n; j++) {
807: if (indexn[j] < 0) { idx++; continue;}
808: #if defined(PETSC_USE_DEBUG)
809: 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);
810: #endif
811: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
812: }
813: }
814: } else {
815: for (i=0; i<m; i++) {
816: if (indexm[i] < 0) { idx += n; continue;}
817: #if defined(PETSC_USE_DEBUG)
818: 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);
819: #endif
820: for (j=0; j<n; j++) {
821: if (indexn[j] < 0) { idx++; continue;}
822: #if defined(PETSC_USE_DEBUG)
823: 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);
824: #endif
825: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
826: }
827: }
828: }
829: }
830: return(0);
831: }
835: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
836: {
837: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
838: PetscInt i,j;
841: /* row-oriented output */
842: for (i=0; i<m; i++) {
843: if (indexm[i] < 0) {v += n;continue;}
844: 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);
845: for (j=0; j<n; j++) {
846: if (indexn[j] < 0) {v++; continue;}
847: 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);
848: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
849: }
850: }
851: return(0);
852: }
854: /* -----------------------------------------------------------------*/
858: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
859: {
860: Mat_SeqDense *a;
862: PetscInt *scols,i,j,nz,header[4];
863: int fd;
864: PetscMPIInt size;
865: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
866: PetscScalar *vals,*svals,*v,*w;
867: MPI_Comm comm;
870: /* force binary viewer to load .info file if it has not yet done so */
871: PetscViewerSetUp(viewer);
872: PetscObjectGetComm((PetscObject)viewer,&comm);
873: MPI_Comm_size(comm,&size);
874: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
875: PetscViewerBinaryGetDescriptor(viewer,&fd);
876: PetscBinaryRead(fd,header,4,PETSC_INT);
877: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
878: M = header[1]; N = header[2]; nz = header[3];
880: /* set global size if not set already*/
881: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
882: MatSetSizes(newmat,M,N,M,N);
883: } else {
884: /* if sizes and type are already set, check if the vector global sizes are correct */
885: MatGetSize(newmat,&grows,&gcols);
886: 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);
887: }
888: a = (Mat_SeqDense*)newmat->data;
889: if (!a->user_alloc) {
890: MatSeqDenseSetPreallocation(newmat,NULL);
891: }
893: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
894: a = (Mat_SeqDense*)newmat->data;
895: v = a->v;
896: /* Allocate some temp space to read in the values and then flip them
897: from row major to column major */
898: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
899: /* read in nonzero values */
900: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
901: /* now flip the values and store them in the matrix*/
902: for (j=0; j<N; j++) {
903: for (i=0; i<M; i++) {
904: *v++ =w[i*N+j];
905: }
906: }
907: PetscFree(w);
908: } else {
909: /* read row lengths */
910: PetscMalloc1(M+1,&rowlengths);
911: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
913: a = (Mat_SeqDense*)newmat->data;
914: v = a->v;
916: /* read column indices and nonzeros */
917: PetscMalloc1(nz+1,&scols);
918: cols = scols;
919: PetscBinaryRead(fd,cols,nz,PETSC_INT);
920: PetscMalloc1(nz+1,&svals);
921: vals = svals;
922: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
924: /* insert into matrix */
925: for (i=0; i<M; i++) {
926: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
927: svals += rowlengths[i]; scols += rowlengths[i];
928: }
929: PetscFree(vals);
930: PetscFree(cols);
931: PetscFree(rowlengths);
932: }
933: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
934: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
935: return(0);
936: }
940: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
941: {
942: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
943: PetscErrorCode ierr;
944: PetscInt i,j;
945: const char *name;
946: PetscScalar *v;
947: PetscViewerFormat format;
948: #if defined(PETSC_USE_COMPLEX)
949: PetscBool allreal = PETSC_TRUE;
950: #endif
953: PetscViewerGetFormat(viewer,&format);
954: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
955: return(0); /* do nothing for now */
956: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
957: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
958: for (i=0; i<A->rmap->n; i++) {
959: v = a->v + i;
960: PetscViewerASCIIPrintf(viewer,"row %D:",i);
961: for (j=0; j<A->cmap->n; j++) {
962: #if defined(PETSC_USE_COMPLEX)
963: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
964: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
965: } else if (PetscRealPart(*v)) {
966: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
967: }
968: #else
969: if (*v) {
970: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
971: }
972: #endif
973: v += a->lda;
974: }
975: PetscViewerASCIIPrintf(viewer,"\n");
976: }
977: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
978: } else {
979: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
980: #if defined(PETSC_USE_COMPLEX)
981: /* determine if matrix has all real values */
982: v = a->v;
983: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
984: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
985: }
986: #endif
987: if (format == PETSC_VIEWER_ASCII_MATLAB) {
988: PetscObjectGetName((PetscObject)A,&name);
989: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
990: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
991: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
992: }
994: for (i=0; i<A->rmap->n; i++) {
995: v = a->v + i;
996: for (j=0; j<A->cmap->n; j++) {
997: #if defined(PETSC_USE_COMPLEX)
998: if (allreal) {
999: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1000: } else {
1001: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1002: }
1003: #else
1004: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1005: #endif
1006: v += a->lda;
1007: }
1008: PetscViewerASCIIPrintf(viewer,"\n");
1009: }
1010: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1011: PetscViewerASCIIPrintf(viewer,"];\n");
1012: }
1013: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1014: }
1015: PetscViewerFlush(viewer);
1016: return(0);
1017: }
1021: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1022: {
1023: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1024: PetscErrorCode ierr;
1025: int fd;
1026: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1027: PetscScalar *v,*anonz,*vals;
1028: PetscViewerFormat format;
1031: PetscViewerBinaryGetDescriptor(viewer,&fd);
1033: PetscViewerGetFormat(viewer,&format);
1034: if (format == PETSC_VIEWER_NATIVE) {
1035: /* store the matrix as a dense matrix */
1036: PetscMalloc1(4,&col_lens);
1038: col_lens[0] = MAT_FILE_CLASSID;
1039: col_lens[1] = m;
1040: col_lens[2] = n;
1041: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1043: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1044: PetscFree(col_lens);
1046: /* write out matrix, by rows */
1047: PetscMalloc1(m*n+1,&vals);
1048: v = a->v;
1049: for (j=0; j<n; j++) {
1050: for (i=0; i<m; i++) {
1051: vals[j + i*n] = *v++;
1052: }
1053: }
1054: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1055: PetscFree(vals);
1056: } else {
1057: PetscMalloc1(4+nz,&col_lens);
1059: col_lens[0] = MAT_FILE_CLASSID;
1060: col_lens[1] = m;
1061: col_lens[2] = n;
1062: col_lens[3] = nz;
1064: /* store lengths of each row and write (including header) to file */
1065: for (i=0; i<m; i++) col_lens[4+i] = n;
1066: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1068: /* Possibly should write in smaller increments, not whole matrix at once? */
1069: /* store column indices (zero start index) */
1070: ict = 0;
1071: for (i=0; i<m; i++) {
1072: for (j=0; j<n; j++) col_lens[ict++] = j;
1073: }
1074: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1075: PetscFree(col_lens);
1077: /* store nonzero values */
1078: PetscMalloc1(nz+1,&anonz);
1079: ict = 0;
1080: for (i=0; i<m; i++) {
1081: v = a->v + i;
1082: for (j=0; j<n; j++) {
1083: anonz[ict++] = *v; v += a->lda;
1084: }
1085: }
1086: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1087: PetscFree(anonz);
1088: }
1089: return(0);
1090: }
1092: #include <petscdraw.h>
1095: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1096: {
1097: Mat A = (Mat) Aa;
1098: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1099: PetscErrorCode ierr;
1100: PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j;
1101: PetscScalar *v = a->v;
1102: PetscViewer viewer;
1103: PetscDraw popup;
1104: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1105: PetscViewerFormat format;
1108: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1109: PetscViewerGetFormat(viewer,&format);
1110: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1112: /* Loop over matrix elements drawing boxes */
1113: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1114: /* Blue for negative and Red for positive */
1115: color = PETSC_DRAW_BLUE;
1116: for (j = 0; j < n; j++) {
1117: x_l = j;
1118: x_r = x_l + 1.0;
1119: for (i = 0; i < m; i++) {
1120: y_l = m - i - 1.0;
1121: y_r = y_l + 1.0;
1122: if (PetscRealPart(v[j*m+i]) > 0.) {
1123: color = PETSC_DRAW_RED;
1124: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1125: color = PETSC_DRAW_BLUE;
1126: } else {
1127: continue;
1128: }
1129: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1130: }
1131: }
1132: } else {
1133: /* use contour shading to indicate magnitude of values */
1134: /* first determine max of all nonzero values */
1135: for (i = 0; i < m*n; i++) {
1136: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1137: }
1138: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1139: PetscDrawGetPopup(draw,&popup);
1140: if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1141: for (j = 0; j < n; j++) {
1142: x_l = j;
1143: x_r = x_l + 1.0;
1144: for (i = 0; i < m; i++) {
1145: y_l = m - i - 1.0;
1146: y_r = y_l + 1.0;
1147: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1148: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1149: }
1150: }
1151: }
1152: return(0);
1153: }
1157: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1158: {
1159: PetscDraw draw;
1160: PetscBool isnull;
1161: PetscReal xr,yr,xl,yl,h,w;
1165: PetscViewerDrawGetDraw(viewer,0,&draw);
1166: PetscDrawIsNull(draw,&isnull);
1167: if (isnull) return(0);
1169: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1170: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1171: xr += w; yr += h; xl = -w; yl = -h;
1172: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1173: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1174: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1175: return(0);
1176: }
1180: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1181: {
1183: PetscBool iascii,isbinary,isdraw;
1186: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1187: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1188: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1190: if (iascii) {
1191: MatView_SeqDense_ASCII(A,viewer);
1192: } else if (isbinary) {
1193: MatView_SeqDense_Binary(A,viewer);
1194: } else if (isdraw) {
1195: MatView_SeqDense_Draw(A,viewer);
1196: }
1197: return(0);
1198: }
1202: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1203: {
1204: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1208: #if defined(PETSC_USE_LOG)
1209: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1210: #endif
1211: PetscFree(l->pivots);
1212: if (!l->user_alloc) {PetscFree(l->v);}
1213: PetscFree(mat->data);
1215: PetscObjectChangeTypeName((PetscObject)mat,0);
1216: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1217: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1218: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1219: #if defined(PETSC_HAVE_ELEMENTAL)
1220: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1221: #endif
1222: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1223: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1224: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1225: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1226: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1227: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1228: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1229: return(0);
1230: }
1234: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1235: {
1236: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1238: PetscInt k,j,m,n,M;
1239: PetscScalar *v,tmp;
1242: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1243: if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1244: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1245: else {
1246: for (j=0; j<m; j++) {
1247: for (k=0; k<j; k++) {
1248: tmp = v[j + k*M];
1249: v[j + k*M] = v[k + j*M];
1250: v[k + j*M] = tmp;
1251: }
1252: }
1253: }
1254: } else { /* out-of-place transpose */
1255: Mat tmat;
1256: Mat_SeqDense *tmatd;
1257: PetscScalar *v2;
1258: PetscInt M2;
1260: if (reuse == MAT_INITIAL_MATRIX) {
1261: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1262: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1263: MatSetType(tmat,((PetscObject)A)->type_name);
1264: MatSeqDenseSetPreallocation(tmat,NULL);
1265: } else {
1266: tmat = *matout;
1267: }
1268: tmatd = (Mat_SeqDense*)tmat->data;
1269: v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1270: for (j=0; j<n; j++) {
1271: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1272: }
1273: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1274: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1275: *matout = tmat;
1276: }
1277: return(0);
1278: }
1282: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1283: {
1284: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1285: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1286: PetscInt i,j;
1287: PetscScalar *v1,*v2;
1290: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1291: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1292: for (i=0; i<A1->rmap->n; i++) {
1293: v1 = mat1->v+i; v2 = mat2->v+i;
1294: for (j=0; j<A1->cmap->n; j++) {
1295: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1296: v1 += mat1->lda; v2 += mat2->lda;
1297: }
1298: }
1299: *flg = PETSC_TRUE;
1300: return(0);
1301: }
1305: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1306: {
1307: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1309: PetscInt i,n,len;
1310: PetscScalar *x,zero = 0.0;
1313: VecSet(v,zero);
1314: VecGetSize(v,&n);
1315: VecGetArray(v,&x);
1316: len = PetscMin(A->rmap->n,A->cmap->n);
1317: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1318: for (i=0; i<len; i++) {
1319: x[i] = mat->v[i*mat->lda + i];
1320: }
1321: VecRestoreArray(v,&x);
1322: return(0);
1323: }
1327: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1328: {
1329: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1330: const PetscScalar *l,*r;
1331: PetscScalar x,*v;
1332: PetscErrorCode ierr;
1333: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1336: if (ll) {
1337: VecGetSize(ll,&m);
1338: VecGetArrayRead(ll,&l);
1339: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1340: for (i=0; i<m; i++) {
1341: x = l[i];
1342: v = mat->v + i;
1343: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1344: }
1345: VecRestoreArrayRead(ll,&l);
1346: PetscLogFlops(1.0*n*m);
1347: }
1348: if (rr) {
1349: VecGetSize(rr,&n);
1350: VecGetArrayRead(rr,&r);
1351: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1352: for (i=0; i<n; i++) {
1353: x = r[i];
1354: v = mat->v + i*m;
1355: for (j=0; j<m; j++) (*v++) *= x;
1356: }
1357: VecRestoreArrayRead(rr,&r);
1358: PetscLogFlops(1.0*n*m);
1359: }
1360: return(0);
1361: }
1365: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1366: {
1367: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1368: PetscScalar *v = mat->v;
1369: PetscReal sum = 0.0;
1370: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1374: if (type == NORM_FROBENIUS) {
1375: if (lda>m) {
1376: for (j=0; j<A->cmap->n; j++) {
1377: v = mat->v+j*lda;
1378: for (i=0; i<m; i++) {
1379: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1380: }
1381: }
1382: } else {
1383: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1384: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1385: }
1386: }
1387: *nrm = PetscSqrtReal(sum);
1388: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1389: } else if (type == NORM_1) {
1390: *nrm = 0.0;
1391: for (j=0; j<A->cmap->n; j++) {
1392: v = mat->v + j*mat->lda;
1393: sum = 0.0;
1394: for (i=0; i<A->rmap->n; i++) {
1395: sum += PetscAbsScalar(*v); v++;
1396: }
1397: if (sum > *nrm) *nrm = sum;
1398: }
1399: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1400: } else if (type == NORM_INFINITY) {
1401: *nrm = 0.0;
1402: for (j=0; j<A->rmap->n; j++) {
1403: v = mat->v + j;
1404: sum = 0.0;
1405: for (i=0; i<A->cmap->n; i++) {
1406: sum += PetscAbsScalar(*v); v += mat->lda;
1407: }
1408: if (sum > *nrm) *nrm = sum;
1409: }
1410: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1411: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1412: return(0);
1413: }
1417: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1418: {
1419: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1423: switch (op) {
1424: case MAT_ROW_ORIENTED:
1425: aij->roworiented = flg;
1426: break;
1427: case MAT_NEW_NONZERO_LOCATIONS:
1428: case MAT_NEW_NONZERO_LOCATION_ERR:
1429: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1430: case MAT_NEW_DIAGONALS:
1431: case MAT_KEEP_NONZERO_PATTERN:
1432: case MAT_IGNORE_OFF_PROC_ENTRIES:
1433: case MAT_USE_HASH_TABLE:
1434: case MAT_IGNORE_LOWER_TRIANGULAR:
1435: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1436: break;
1437: case MAT_SPD:
1438: case MAT_SYMMETRIC:
1439: case MAT_STRUCTURALLY_SYMMETRIC:
1440: case MAT_HERMITIAN:
1441: case MAT_SYMMETRY_ETERNAL:
1442: /* These options are handled directly by MatSetOption() */
1443: break;
1444: default:
1445: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1446: }
1447: return(0);
1448: }
1452: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1453: {
1454: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1456: PetscInt lda=l->lda,m=A->rmap->n,j;
1459: if (lda>m) {
1460: for (j=0; j<A->cmap->n; j++) {
1461: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1462: }
1463: } else {
1464: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1465: }
1466: return(0);
1467: }
1471: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1472: {
1473: PetscErrorCode ierr;
1474: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1475: PetscInt m = l->lda, n = A->cmap->n, i,j;
1476: PetscScalar *slot,*bb;
1477: const PetscScalar *xx;
1480: #if defined(PETSC_USE_DEBUG)
1481: for (i=0; i<N; i++) {
1482: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1483: 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);
1484: }
1485: #endif
1487: /* fix right hand side if needed */
1488: if (x && b) {
1489: VecGetArrayRead(x,&xx);
1490: VecGetArray(b,&bb);
1491: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1492: VecRestoreArrayRead(x,&xx);
1493: VecRestoreArray(b,&bb);
1494: }
1496: for (i=0; i<N; i++) {
1497: slot = l->v + rows[i];
1498: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1499: }
1500: if (diag != 0.0) {
1501: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1502: for (i=0; i<N; i++) {
1503: slot = l->v + (m+1)*rows[i];
1504: *slot = diag;
1505: }
1506: }
1507: return(0);
1508: }
1512: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1513: {
1514: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1517: 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");
1518: *array = mat->v;
1519: return(0);
1520: }
1524: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1525: {
1527: *array = 0; /* user cannot accidently use the array later */
1528: return(0);
1529: }
1533: /*@C
1534: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1536: Not Collective
1538: Input Parameter:
1539: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1541: Output Parameter:
1542: . array - pointer to the data
1544: Level: intermediate
1546: .seealso: MatDenseRestoreArray()
1547: @*/
1548: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1549: {
1553: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1554: return(0);
1555: }
1559: /*@C
1560: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1562: Not Collective
1564: Input Parameters:
1565: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1566: . array - pointer to the data
1568: Level: intermediate
1570: .seealso: MatDenseGetArray()
1571: @*/
1572: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1573: {
1577: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1578: return(0);
1579: }
1583: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1584: {
1585: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1587: PetscInt i,j,nrows,ncols;
1588: const PetscInt *irow,*icol;
1589: PetscScalar *av,*bv,*v = mat->v;
1590: Mat newmat;
1593: ISGetIndices(isrow,&irow);
1594: ISGetIndices(iscol,&icol);
1595: ISGetLocalSize(isrow,&nrows);
1596: ISGetLocalSize(iscol,&ncols);
1598: /* Check submatrixcall */
1599: if (scall == MAT_REUSE_MATRIX) {
1600: PetscInt n_cols,n_rows;
1601: MatGetSize(*B,&n_rows,&n_cols);
1602: if (n_rows != nrows || n_cols != ncols) {
1603: /* resize the result matrix to match number of requested rows/columns */
1604: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1605: }
1606: newmat = *B;
1607: } else {
1608: /* Create and fill new matrix */
1609: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1610: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1611: MatSetType(newmat,((PetscObject)A)->type_name);
1612: MatSeqDenseSetPreallocation(newmat,NULL);
1613: }
1615: /* Now extract the data pointers and do the copy,column at a time */
1616: bv = ((Mat_SeqDense*)newmat->data)->v;
1618: for (i=0; i<ncols; i++) {
1619: av = v + mat->lda*icol[i];
1620: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1621: }
1623: /* Assemble the matrices so that the correct flags are set */
1624: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1625: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1627: /* Free work space */
1628: ISRestoreIndices(isrow,&irow);
1629: ISRestoreIndices(iscol,&icol);
1630: *B = newmat;
1631: return(0);
1632: }
1636: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1637: {
1639: PetscInt i;
1642: if (scall == MAT_INITIAL_MATRIX) {
1643: PetscMalloc1(n+1,B);
1644: }
1646: for (i=0; i<n; i++) {
1647: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1648: }
1649: return(0);
1650: }
1654: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1655: {
1657: return(0);
1658: }
1662: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1663: {
1665: return(0);
1666: }
1670: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1671: {
1672: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1674: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1677: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1678: if (A->ops->copy != B->ops->copy) {
1679: MatCopy_Basic(A,B,str);
1680: return(0);
1681: }
1682: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1683: if (lda1>m || lda2>m) {
1684: for (j=0; j<n; j++) {
1685: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1686: }
1687: } else {
1688: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1689: }
1690: return(0);
1691: }
1695: PetscErrorCode MatSetUp_SeqDense(Mat A)
1696: {
1700: MatSeqDenseSetPreallocation(A,0);
1701: return(0);
1702: }
1706: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1707: {
1708: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1709: PetscInt i,nz = A->rmap->n*A->cmap->n;
1710: PetscScalar *aa = a->v;
1713: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1714: return(0);
1715: }
1719: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1720: {
1721: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1722: PetscInt i,nz = A->rmap->n*A->cmap->n;
1723: PetscScalar *aa = a->v;
1726: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1727: return(0);
1728: }
1732: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1733: {
1734: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1735: PetscInt i,nz = A->rmap->n*A->cmap->n;
1736: PetscScalar *aa = a->v;
1739: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1740: return(0);
1741: }
1743: /* ----------------------------------------------------------------*/
1746: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1747: {
1751: if (scall == MAT_INITIAL_MATRIX) {
1752: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1753: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1754: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1755: }
1756: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1757: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1758: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1759: return(0);
1760: }
1764: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1765: {
1767: PetscInt m=A->rmap->n,n=B->cmap->n;
1768: Mat Cmat;
1771: 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);
1772: MatCreate(PETSC_COMM_SELF,&Cmat);
1773: MatSetSizes(Cmat,m,n,m,n);
1774: MatSetType(Cmat,MATSEQDENSE);
1775: MatSeqDenseSetPreallocation(Cmat,NULL);
1777: *C = Cmat;
1778: return(0);
1779: }
1783: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1784: {
1785: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1786: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1787: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1788: PetscBLASInt m,n,k;
1789: PetscScalar _DOne=1.0,_DZero=0.0;
1793: PetscBLASIntCast(A->rmap->n,&m);
1794: PetscBLASIntCast(B->cmap->n,&n);
1795: PetscBLASIntCast(A->cmap->n,&k);
1796: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1797: return(0);
1798: }
1802: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1803: {
1807: if (scall == MAT_INITIAL_MATRIX) {
1808: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1809: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1810: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1811: }
1812: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1813: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1814: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1815: return(0);
1816: }
1820: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1821: {
1823: PetscInt m=A->cmap->n,n=B->cmap->n;
1824: Mat Cmat;
1827: 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);
1828: MatCreate(PETSC_COMM_SELF,&Cmat);
1829: MatSetSizes(Cmat,m,n,m,n);
1830: MatSetType(Cmat,MATSEQDENSE);
1831: MatSeqDenseSetPreallocation(Cmat,NULL);
1833: Cmat->assembled = PETSC_TRUE;
1835: *C = Cmat;
1836: return(0);
1837: }
1841: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1842: {
1843: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1844: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1845: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1846: PetscBLASInt m,n,k;
1847: PetscScalar _DOne=1.0,_DZero=0.0;
1851: PetscBLASIntCast(A->cmap->n,&m);
1852: PetscBLASIntCast(B->cmap->n,&n);
1853: PetscBLASIntCast(A->rmap->n,&k);
1854: /*
1855: Note the m and n arguments below are the number rows and columns of A', not A!
1856: */
1857: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1858: return(0);
1859: }
1863: PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1864: {
1865: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1867: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1868: PetscScalar *x;
1869: MatScalar *aa = a->v;
1872: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1874: VecSet(v,0.0);
1875: VecGetArray(v,&x);
1876: VecGetLocalSize(v,&p);
1877: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1878: for (i=0; i<m; i++) {
1879: x[i] = aa[i]; if (idx) idx[i] = 0;
1880: for (j=1; j<n; j++) {
1881: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1882: }
1883: }
1884: VecRestoreArray(v,&x);
1885: return(0);
1886: }
1890: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1891: {
1892: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1894: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1895: PetscScalar *x;
1896: PetscReal atmp;
1897: MatScalar *aa = a->v;
1900: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1902: VecSet(v,0.0);
1903: VecGetArray(v,&x);
1904: VecGetLocalSize(v,&p);
1905: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1906: for (i=0; i<m; i++) {
1907: x[i] = PetscAbsScalar(aa[i]);
1908: for (j=1; j<n; j++) {
1909: atmp = PetscAbsScalar(aa[i+m*j]);
1910: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1911: }
1912: }
1913: VecRestoreArray(v,&x);
1914: return(0);
1915: }
1919: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1920: {
1921: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1923: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1924: PetscScalar *x;
1925: MatScalar *aa = a->v;
1928: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1930: VecSet(v,0.0);
1931: VecGetArray(v,&x);
1932: VecGetLocalSize(v,&p);
1933: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1934: for (i=0; i<m; i++) {
1935: x[i] = aa[i]; if (idx) idx[i] = 0;
1936: for (j=1; j<n; j++) {
1937: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1938: }
1939: }
1940: VecRestoreArray(v,&x);
1941: return(0);
1942: }
1946: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1947: {
1948: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1950: PetscScalar *x;
1953: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1955: VecGetArray(v,&x);
1956: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1957: VecRestoreArray(v,&x);
1958: return(0);
1959: }
1964: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1965: {
1967: PetscInt i,j,m,n;
1968: PetscScalar *a;
1971: MatGetSize(A,&m,&n);
1972: PetscMemzero(norms,n*sizeof(PetscReal));
1973: MatDenseGetArray(A,&a);
1974: if (type == NORM_2) {
1975: for (i=0; i<n; i++) {
1976: for (j=0; j<m; j++) {
1977: norms[i] += PetscAbsScalar(a[j]*a[j]);
1978: }
1979: a += m;
1980: }
1981: } else if (type == NORM_1) {
1982: for (i=0; i<n; i++) {
1983: for (j=0; j<m; j++) {
1984: norms[i] += PetscAbsScalar(a[j]);
1985: }
1986: a += m;
1987: }
1988: } else if (type == NORM_INFINITY) {
1989: for (i=0; i<n; i++) {
1990: for (j=0; j<m; j++) {
1991: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1992: }
1993: a += m;
1994: }
1995: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1996: MatDenseRestoreArray(A,&a);
1997: if (type == NORM_2) {
1998: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1999: }
2000: return(0);
2001: }
2005: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2006: {
2008: PetscScalar *a;
2009: PetscInt m,n,i;
2012: MatGetSize(x,&m,&n);
2013: MatDenseGetArray(x,&a);
2014: for (i=0; i<m*n; i++) {
2015: PetscRandomGetValue(rctx,a+i);
2016: }
2017: MatDenseRestoreArray(x,&a);
2018: return(0);
2019: }
2022: /* -------------------------------------------------------------------*/
2023: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2024: MatGetRow_SeqDense,
2025: MatRestoreRow_SeqDense,
2026: MatMult_SeqDense,
2027: /* 4*/ MatMultAdd_SeqDense,
2028: MatMultTranspose_SeqDense,
2029: MatMultTransposeAdd_SeqDense,
2030: 0,
2031: 0,
2032: 0,
2033: /* 10*/ 0,
2034: MatLUFactor_SeqDense,
2035: MatCholeskyFactor_SeqDense,
2036: MatSOR_SeqDense,
2037: MatTranspose_SeqDense,
2038: /* 15*/ MatGetInfo_SeqDense,
2039: MatEqual_SeqDense,
2040: MatGetDiagonal_SeqDense,
2041: MatDiagonalScale_SeqDense,
2042: MatNorm_SeqDense,
2043: /* 20*/ MatAssemblyBegin_SeqDense,
2044: MatAssemblyEnd_SeqDense,
2045: MatSetOption_SeqDense,
2046: MatZeroEntries_SeqDense,
2047: /* 24*/ MatZeroRows_SeqDense,
2048: 0,
2049: 0,
2050: 0,
2051: 0,
2052: /* 29*/ MatSetUp_SeqDense,
2053: 0,
2054: 0,
2055: 0,
2056: 0,
2057: /* 34*/ MatDuplicate_SeqDense,
2058: 0,
2059: 0,
2060: 0,
2061: 0,
2062: /* 39*/ MatAXPY_SeqDense,
2063: MatGetSubMatrices_SeqDense,
2064: 0,
2065: MatGetValues_SeqDense,
2066: MatCopy_SeqDense,
2067: /* 44*/ MatGetRowMax_SeqDense,
2068: MatScale_SeqDense,
2069: MatShift_Basic,
2070: 0,
2071: 0,
2072: /* 49*/ MatSetRandom_SeqDense,
2073: 0,
2074: 0,
2075: 0,
2076: 0,
2077: /* 54*/ 0,
2078: 0,
2079: 0,
2080: 0,
2081: 0,
2082: /* 59*/ 0,
2083: MatDestroy_SeqDense,
2084: MatView_SeqDense,
2085: 0,
2086: 0,
2087: /* 64*/ 0,
2088: 0,
2089: 0,
2090: 0,
2091: 0,
2092: /* 69*/ MatGetRowMaxAbs_SeqDense,
2093: 0,
2094: 0,
2095: 0,
2096: 0,
2097: /* 74*/ 0,
2098: 0,
2099: 0,
2100: 0,
2101: 0,
2102: /* 79*/ 0,
2103: 0,
2104: 0,
2105: 0,
2106: /* 83*/ MatLoad_SeqDense,
2107: 0,
2108: MatIsHermitian_SeqDense,
2109: 0,
2110: 0,
2111: 0,
2112: /* 89*/ MatMatMult_SeqDense_SeqDense,
2113: MatMatMultSymbolic_SeqDense_SeqDense,
2114: MatMatMultNumeric_SeqDense_SeqDense,
2115: 0,
2116: 0,
2117: /* 94*/ 0,
2118: 0,
2119: 0,
2120: 0,
2121: 0,
2122: /* 99*/ 0,
2123: 0,
2124: 0,
2125: MatConjugate_SeqDense,
2126: 0,
2127: /*104*/ 0,
2128: MatRealPart_SeqDense,
2129: MatImaginaryPart_SeqDense,
2130: 0,
2131: 0,
2132: /*109*/ MatMatSolve_SeqDense,
2133: 0,
2134: MatGetRowMin_SeqDense,
2135: MatGetColumnVector_SeqDense,
2136: 0,
2137: /*114*/ 0,
2138: 0,
2139: 0,
2140: 0,
2141: 0,
2142: /*119*/ 0,
2143: 0,
2144: 0,
2145: 0,
2146: 0,
2147: /*124*/ 0,
2148: MatGetColumnNorms_SeqDense,
2149: 0,
2150: 0,
2151: 0,
2152: /*129*/ 0,
2153: MatTransposeMatMult_SeqDense_SeqDense,
2154: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2155: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2156: 0,
2157: /*134*/ 0,
2158: 0,
2159: 0,
2160: 0,
2161: 0,
2162: /*139*/ 0,
2163: 0,
2164: 0
2165: };
2169: /*@C
2170: MatCreateSeqDense - Creates a sequential dense matrix that
2171: is stored in column major order (the usual Fortran 77 manner). Many
2172: of the matrix operations use the BLAS and LAPACK routines.
2174: Collective on MPI_Comm
2176: Input Parameters:
2177: + comm - MPI communicator, set to PETSC_COMM_SELF
2178: . m - number of rows
2179: . n - number of columns
2180: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2181: to control all matrix memory allocation.
2183: Output Parameter:
2184: . A - the matrix
2186: Notes:
2187: The data input variable is intended primarily for Fortran programmers
2188: who wish to allocate their own matrix memory space. Most users should
2189: set data=NULL.
2191: Level: intermediate
2193: .keywords: dense, matrix, LAPACK, BLAS
2195: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2196: @*/
2197: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2198: {
2202: MatCreate(comm,A);
2203: MatSetSizes(*A,m,n,m,n);
2204: MatSetType(*A,MATSEQDENSE);
2205: MatSeqDenseSetPreallocation(*A,data);
2206: return(0);
2207: }
2211: /*@C
2212: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2214: Collective on MPI_Comm
2216: Input Parameters:
2217: + B - the matrix
2218: - data - the array (or NULL)
2220: Notes:
2221: The data input variable is intended primarily for Fortran programmers
2222: who wish to allocate their own matrix memory space. Most users should
2223: need not call this routine.
2225: Level: intermediate
2227: .keywords: dense, matrix, LAPACK, BLAS
2229: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2231: @*/
2232: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2233: {
2237: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2238: return(0);
2239: }
2243: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2244: {
2245: Mat_SeqDense *b;
2249: B->preallocated = PETSC_TRUE;
2251: PetscLayoutSetUp(B->rmap);
2252: PetscLayoutSetUp(B->cmap);
2254: b = (Mat_SeqDense*)B->data;
2255: b->Mmax = B->rmap->n;
2256: b->Nmax = B->cmap->n;
2257: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2259: if (!data) { /* petsc-allocated storage */
2260: if (!b->user_alloc) { PetscFree(b->v); }
2261: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2262: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2264: b->user_alloc = PETSC_FALSE;
2265: } else { /* user-allocated storage */
2266: if (!b->user_alloc) { PetscFree(b->v); }
2267: b->v = data;
2268: b->user_alloc = PETSC_TRUE;
2269: }
2270: B->assembled = PETSC_TRUE;
2271: return(0);
2272: }
2274: #if defined(PETSC_HAVE_ELEMENTAL)
2277: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2278: {
2279: Mat mat_elemental;
2281: PetscScalar *array,*v_colwise;
2282: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2285: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2286: MatDenseGetArray(A,&array);
2287: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2288: k = 0;
2289: for (j=0; j<N; j++) {
2290: cols[j] = j;
2291: for (i=0; i<M; i++) {
2292: v_colwise[j*M+i] = array[k++];
2293: }
2294: }
2295: for (i=0; i<M; i++) {
2296: rows[i] = i;
2297: }
2298: MatDenseRestoreArray(A,&array);
2300: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2301: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2302: MatSetType(mat_elemental,MATELEMENTAL);
2303: MatSetUp(mat_elemental);
2305: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2306: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2307: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2308: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2309: PetscFree3(v_colwise,rows,cols);
2311: if (reuse == MAT_REUSE_MATRIX) {
2312: MatHeaderReplace(A,mat_elemental);
2313: } else {
2314: *newmat = mat_elemental;
2315: }
2316: return(0);
2317: }
2318: #endif
2322: /*@C
2323: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2325: Input parameter:
2326: + A - the matrix
2327: - lda - the leading dimension
2329: Notes:
2330: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2331: it asserts that the preallocation has a leading dimension (the LDA parameter
2332: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2334: Level: intermediate
2336: .keywords: dense, matrix, LAPACK, BLAS
2338: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2340: @*/
2341: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2342: {
2343: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2346: 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);
2347: b->lda = lda;
2348: b->changelda = PETSC_FALSE;
2349: b->Mmax = PetscMax(b->Mmax,lda);
2350: return(0);
2351: }
2353: /*MC
2354: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2356: Options Database Keys:
2357: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2359: Level: beginner
2361: .seealso: MatCreateSeqDense()
2363: M*/
2367: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2368: {
2369: Mat_SeqDense *b;
2371: PetscMPIInt size;
2374: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2375: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2377: PetscNewLog(B,&b);
2378: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2379: B->data = (void*)b;
2381: b->pivots = 0;
2382: b->roworiented = PETSC_TRUE;
2383: b->v = 0;
2384: b->changelda = PETSC_FALSE;
2386: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2387: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2388: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2389: #if defined(PETSC_HAVE_ELEMENTAL)
2390: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2391: #endif
2392: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2393: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2394: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2395: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2396: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2397: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2398: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2399: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2400: return(0);
2401: }