Actual source code: dense.c
petsc-3.5.4 2015-05-23
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_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
14: {
15: Mat B;
16: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
18: PetscInt i, j;
19: PetscInt *rows, *nnz;
20: MatScalar *aa = a->v, *vals;
23: MatCreate(PetscObjectComm((PetscObject)A),&B);
24: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
25: MatSetType(B,MATSEQAIJ);
26: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
27: for (j=0; j<A->cmap->n; j++) {
28: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
29: aa += a->lda;
30: }
31: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
32: aa = a->v;
33: for (j=0; j<A->cmap->n; j++) {
34: PetscInt numRows = 0;
35: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
36: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
37: aa += a->lda;
38: }
39: PetscFree3(rows,nnz,vals);
40: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
41: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
43: if (reuse == MAT_REUSE_MATRIX) {
44: MatHeaderReplace(A,B);
45: } else {
46: *newmat = B;
47: }
48: return(0);
49: }
53: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
54: {
55: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
56: PetscScalar oalpha = alpha;
57: PetscInt j;
58: PetscBLASInt N,m,ldax,lday,one = 1;
62: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
63: PetscBLASIntCast(X->rmap->n,&m);
64: PetscBLASIntCast(x->lda,&ldax);
65: PetscBLASIntCast(y->lda,&lday);
66: if (ldax>m || lday>m) {
67: for (j=0; j<X->cmap->n; j++) {
68: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
69: }
70: } else {
71: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
72: }
73: PetscObjectStateIncrease((PetscObject)Y);
74: PetscLogFlops(PetscMax(2*N-1,0));
75: return(0);
76: }
80: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
81: {
82: PetscInt N = A->rmap->n*A->cmap->n;
85: info->block_size = 1.0;
86: info->nz_allocated = (double)N;
87: info->nz_used = (double)N;
88: info->nz_unneeded = (double)0;
89: info->assemblies = (double)A->num_ass;
90: info->mallocs = 0;
91: info->memory = ((PetscObject)A)->mem;
92: info->fill_ratio_given = 0;
93: info->fill_ratio_needed = 0;
94: info->factor_mallocs = 0;
95: return(0);
96: }
100: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
101: {
102: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
103: PetscScalar oalpha = alpha;
105: PetscBLASInt one = 1,j,nz,lda;
108: PetscBLASIntCast(a->lda,&lda);
109: if (lda>A->rmap->n) {
110: PetscBLASIntCast(A->rmap->n,&nz);
111: for (j=0; j<A->cmap->n; j++) {
112: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
113: }
114: } else {
115: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
116: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
117: }
118: PetscLogFlops(nz);
119: return(0);
120: }
124: PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
125: {
126: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
127: PetscInt i,j,m = A->rmap->n,N;
128: PetscScalar *v = a->v;
131: *fl = PETSC_FALSE;
132: if (A->rmap->n != A->cmap->n) return(0);
133: N = a->lda;
135: for (i=0; i<m; i++) {
136: for (j=i+1; j<m; j++) {
137: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
138: }
139: }
140: *fl = PETSC_TRUE;
141: return(0);
142: }
146: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
147: {
148: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
150: PetscInt lda = (PetscInt)mat->lda,j,m;
153: PetscLayoutReference(A->rmap,&newi->rmap);
154: PetscLayoutReference(A->cmap,&newi->cmap);
155: MatSeqDenseSetPreallocation(newi,NULL);
156: if (cpvalues == MAT_COPY_VALUES) {
157: l = (Mat_SeqDense*)newi->data;
158: if (lda>A->rmap->n) {
159: m = A->rmap->n;
160: for (j=0; j<A->cmap->n; j++) {
161: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
162: }
163: } else {
164: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
165: }
166: }
167: newi->assembled = PETSC_TRUE;
168: return(0);
169: }
173: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
174: {
178: MatCreate(PetscObjectComm((PetscObject)A),newmat);
179: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
180: MatSetType(*newmat,((PetscObject)A)->type_name);
181: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
182: return(0);
183: }
186: extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
190: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
191: {
192: MatFactorInfo info;
196: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
197: MatLUFactor_SeqDense(fact,0,0,&info);
198: return(0);
199: }
203: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
204: {
205: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
207: PetscScalar *x,*y;
208: PetscBLASInt one = 1,info,m;
211: PetscBLASIntCast(A->rmap->n,&m);
212: VecGetArray(xx,&x);
213: VecGetArray(yy,&y);
214: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
215: if (A->factortype == MAT_FACTOR_LU) {
216: #if defined(PETSC_MISSING_LAPACK_GETRS)
217: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
218: #else
219: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
220: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
221: #endif
222: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
223: #if defined(PETSC_MISSING_LAPACK_POTRS)
224: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
225: #else
226: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
227: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
228: #endif
229: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
230: VecRestoreArray(xx,&x);
231: VecRestoreArray(yy,&y);
232: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
233: return(0);
234: }
238: PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
239: {
240: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
242: PetscScalar *b,*x;
243: PetscInt n;
244: PetscBLASInt nrhs,info,m;
245: PetscBool flg;
248: PetscBLASIntCast(A->rmap->n,&m);
249: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
250: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
251: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
252: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
254: MatGetSize(B,NULL,&n);
255: PetscBLASIntCast(n,&nrhs);
256: MatDenseGetArray(B,&b);
257: MatDenseGetArray(X,&x);
259: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
261: if (A->factortype == MAT_FACTOR_LU) {
262: #if defined(PETSC_MISSING_LAPACK_GETRS)
263: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
264: #else
265: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
266: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
267: #endif
268: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
269: #if defined(PETSC_MISSING_LAPACK_POTRS)
270: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
271: #else
272: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
273: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
274: #endif
275: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
277: MatDenseRestoreArray(B,&b);
278: MatDenseRestoreArray(X,&x);
279: PetscLogFlops(nrhs*(2.0*m*m - m));
280: return(0);
281: }
285: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
286: {
287: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
289: PetscScalar *x,*y;
290: PetscBLASInt one = 1,info,m;
293: PetscBLASIntCast(A->rmap->n,&m);
294: VecGetArray(xx,&x);
295: VecGetArray(yy,&y);
296: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
297: /* assume if pivots exist then use LU; else Cholesky */
298: if (mat->pivots) {
299: #if defined(PETSC_MISSING_LAPACK_GETRS)
300: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
301: #else
302: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
303: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
304: #endif
305: } else {
306: #if defined(PETSC_MISSING_LAPACK_POTRS)
307: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
308: #else
309: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
310: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
311: #endif
312: }
313: VecRestoreArray(xx,&x);
314: VecRestoreArray(yy,&y);
315: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
316: return(0);
317: }
321: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
322: {
323: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
325: PetscScalar *x,*y,sone = 1.0;
326: Vec tmp = 0;
327: PetscBLASInt one = 1,info,m;
330: PetscBLASIntCast(A->rmap->n,&m);
331: VecGetArray(xx,&x);
332: VecGetArray(yy,&y);
333: if (!A->rmap->n || !A->cmap->n) return(0);
334: if (yy == zz) {
335: VecDuplicate(yy,&tmp);
336: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
337: VecCopy(yy,tmp);
338: }
339: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
340: /* assume if pivots exist then use LU; else Cholesky */
341: if (mat->pivots) {
342: #if defined(PETSC_MISSING_LAPACK_GETRS)
343: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
344: #else
345: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
346: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
347: #endif
348: } else {
349: #if defined(PETSC_MISSING_LAPACK_POTRS)
350: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
351: #else
352: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
353: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
354: #endif
355: }
356: if (tmp) {
357: VecAXPY(yy,sone,tmp);
358: VecDestroy(&tmp);
359: } else {
360: VecAXPY(yy,sone,zz);
361: }
362: VecRestoreArray(xx,&x);
363: VecRestoreArray(yy,&y);
364: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
365: return(0);
366: }
370: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
371: {
372: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
374: PetscScalar *x,*y,sone = 1.0;
375: Vec tmp;
376: PetscBLASInt one = 1,info,m;
379: PetscBLASIntCast(A->rmap->n,&m);
380: if (!A->rmap->n || !A->cmap->n) return(0);
381: VecGetArray(xx,&x);
382: VecGetArray(yy,&y);
383: if (yy == zz) {
384: VecDuplicate(yy,&tmp);
385: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
386: VecCopy(yy,tmp);
387: }
388: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
389: /* assume if pivots exist then use LU; else Cholesky */
390: if (mat->pivots) {
391: #if defined(PETSC_MISSING_LAPACK_GETRS)
392: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
393: #else
394: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
395: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
396: #endif
397: } else {
398: #if defined(PETSC_MISSING_LAPACK_POTRS)
399: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
400: #else
401: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
402: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
403: #endif
404: }
405: if (tmp) {
406: VecAXPY(yy,sone,tmp);
407: VecDestroy(&tmp);
408: } else {
409: VecAXPY(yy,sone,zz);
410: }
411: VecRestoreArray(xx,&x);
412: VecRestoreArray(yy,&y);
413: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
414: return(0);
415: }
417: /* ---------------------------------------------------------------*/
418: /* COMMENT: I have chosen to hide row permutation in the pivots,
419: rather than put it in the Mat->row slot.*/
422: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
423: {
424: #if defined(PETSC_MISSING_LAPACK_GETRF)
426: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
427: #else
428: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
430: PetscBLASInt n,m,info;
433: PetscBLASIntCast(A->cmap->n,&n);
434: PetscBLASIntCast(A->rmap->n,&m);
435: if (!mat->pivots) {
436: PetscMalloc1((A->rmap->n+1),&mat->pivots);
437: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
438: }
439: if (!A->rmap->n || !A->cmap->n) return(0);
440: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
441: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
442: PetscFPTrapPop();
444: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
445: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
446: A->ops->solve = MatSolve_SeqDense;
447: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
448: A->ops->solveadd = MatSolveAdd_SeqDense;
449: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
450: A->factortype = MAT_FACTOR_LU;
452: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
453: #endif
454: return(0);
455: }
459: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
460: {
461: #if defined(PETSC_MISSING_LAPACK_POTRF)
463: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
464: #else
465: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
467: PetscBLASInt info,n;
470: PetscBLASIntCast(A->cmap->n,&n);
471: PetscFree(mat->pivots);
473: if (!A->rmap->n || !A->cmap->n) return(0);
474: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
475: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
476: A->ops->solve = MatSolve_SeqDense;
477: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
478: A->ops->solveadd = MatSolveAdd_SeqDense;
479: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
480: A->factortype = MAT_FACTOR_CHOLESKY;
482: PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
483: #endif
484: return(0);
485: }
490: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
491: {
493: MatFactorInfo info;
496: info.fill = 1.0;
498: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
499: MatCholeskyFactor_SeqDense(fact,0,&info);
500: return(0);
501: }
505: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
506: {
508: fact->assembled = PETSC_TRUE;
509: fact->preallocated = PETSC_TRUE;
510: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
511: return(0);
512: }
516: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
517: {
519: fact->preallocated = PETSC_TRUE;
520: fact->assembled = PETSC_TRUE;
521: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
522: return(0);
523: }
527: PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
528: {
532: MatCreate(PetscObjectComm((PetscObject)A),fact);
533: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
534: MatSetType(*fact,((PetscObject)A)->type_name);
535: if (ftype == MAT_FACTOR_LU) {
536: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
537: } else {
538: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
539: }
540: (*fact)->factortype = ftype;
541: return(0);
542: }
544: /* ------------------------------------------------------------------*/
547: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
548: {
549: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
550: PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt;
552: PetscInt m = A->rmap->n,i;
553: PetscBLASInt o = 1,bm;
556: PetscBLASIntCast(m,&bm);
557: if (flag & SOR_ZERO_INITIAL_GUESS) {
558: /* this is a hack fix, should have another version without the second BLASdot */
559: VecSet(xx,zero);
560: }
561: VecGetArray(xx,&x);
562: VecGetArray(bb,&b);
563: its = its*lits;
564: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
565: while (its--) {
566: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
567: for (i=0; i<m; i++) {
568: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
569: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
570: }
571: }
572: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
573: for (i=m-1; i>=0; i--) {
574: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
575: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
576: }
577: }
578: }
579: VecRestoreArray(bb,&b);
580: VecRestoreArray(xx,&x);
581: return(0);
582: }
584: /* -----------------------------------------------------------------*/
587: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
588: {
589: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
590: PetscScalar *v = mat->v,*x,*y;
592: PetscBLASInt m, n,_One=1;
593: PetscScalar _DOne=1.0,_DZero=0.0;
596: PetscBLASIntCast(A->rmap->n,&m);
597: PetscBLASIntCast(A->cmap->n,&n);
598: if (!A->rmap->n || !A->cmap->n) return(0);
599: VecGetArray(xx,&x);
600: VecGetArray(yy,&y);
601: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
602: VecRestoreArray(xx,&x);
603: VecRestoreArray(yy,&y);
604: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
605: return(0);
606: }
610: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
611: {
612: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
613: PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
615: PetscBLASInt m, n, _One=1;
618: PetscBLASIntCast(A->rmap->n,&m);
619: PetscBLASIntCast(A->cmap->n,&n);
620: if (!A->rmap->n || !A->cmap->n) return(0);
621: VecGetArray(xx,&x);
622: VecGetArray(yy,&y);
623: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
624: VecRestoreArray(xx,&x);
625: VecRestoreArray(yy,&y);
626: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
627: return(0);
628: }
632: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
633: {
634: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
635: PetscScalar *v = mat->v,*x,*y,_DOne=1.0;
637: PetscBLASInt m, n, _One=1;
640: PetscBLASIntCast(A->rmap->n,&m);
641: PetscBLASIntCast(A->cmap->n,&n);
642: if (!A->rmap->n || !A->cmap->n) return(0);
643: if (zz != yy) {VecCopy(zz,yy);}
644: VecGetArray(xx,&x);
645: VecGetArray(yy,&y);
646: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
647: VecRestoreArray(xx,&x);
648: VecRestoreArray(yy,&y);
649: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
650: return(0);
651: }
655: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
656: {
657: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
658: PetscScalar *v = mat->v,*x,*y;
660: PetscBLASInt m, n, _One=1;
661: PetscScalar _DOne=1.0;
664: PetscBLASIntCast(A->rmap->n,&m);
665: PetscBLASIntCast(A->cmap->n,&n);
666: if (!A->rmap->n || !A->cmap->n) return(0);
667: if (zz != yy) {VecCopy(zz,yy);}
668: VecGetArray(xx,&x);
669: VecGetArray(yy,&y);
670: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
671: VecRestoreArray(xx,&x);
672: VecRestoreArray(yy,&y);
673: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
674: return(0);
675: }
677: /* -----------------------------------------------------------------*/
680: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
681: {
682: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
683: PetscScalar *v;
685: PetscInt i;
688: *ncols = A->cmap->n;
689: if (cols) {
690: PetscMalloc1((A->cmap->n+1),cols);
691: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
692: }
693: if (vals) {
694: PetscMalloc1((A->cmap->n+1),vals);
695: v = mat->v + row;
696: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
697: }
698: return(0);
699: }
703: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
704: {
708: if (cols) {PetscFree(*cols);}
709: if (vals) {PetscFree(*vals); }
710: return(0);
711: }
712: /* ----------------------------------------------------------------*/
715: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
716: {
717: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
718: PetscInt i,j,idx=0;
721: if (!mat->roworiented) {
722: if (addv == INSERT_VALUES) {
723: for (j=0; j<n; j++) {
724: if (indexn[j] < 0) {idx += m; continue;}
725: #if defined(PETSC_USE_DEBUG)
726: 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);
727: #endif
728: for (i=0; i<m; i++) {
729: if (indexm[i] < 0) {idx++; continue;}
730: #if defined(PETSC_USE_DEBUG)
731: 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);
732: #endif
733: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
734: }
735: }
736: } else {
737: for (j=0; j<n; j++) {
738: if (indexn[j] < 0) {idx += m; continue;}
739: #if defined(PETSC_USE_DEBUG)
740: 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);
741: #endif
742: for (i=0; i<m; i++) {
743: if (indexm[i] < 0) {idx++; continue;}
744: #if defined(PETSC_USE_DEBUG)
745: 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);
746: #endif
747: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
748: }
749: }
750: }
751: } else {
752: if (addv == INSERT_VALUES) {
753: for (i=0; i<m; i++) {
754: if (indexm[i] < 0) { idx += n; continue;}
755: #if defined(PETSC_USE_DEBUG)
756: 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);
757: #endif
758: for (j=0; j<n; j++) {
759: if (indexn[j] < 0) { idx++; continue;}
760: #if defined(PETSC_USE_DEBUG)
761: 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);
762: #endif
763: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
764: }
765: }
766: } else {
767: for (i=0; i<m; i++) {
768: if (indexm[i] < 0) { idx += n; continue;}
769: #if defined(PETSC_USE_DEBUG)
770: 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);
771: #endif
772: for (j=0; j<n; j++) {
773: if (indexn[j] < 0) { idx++; continue;}
774: #if defined(PETSC_USE_DEBUG)
775: 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);
776: #endif
777: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
778: }
779: }
780: }
781: }
782: return(0);
783: }
787: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
788: {
789: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
790: PetscInt i,j;
793: /* row-oriented output */
794: for (i=0; i<m; i++) {
795: if (indexm[i] < 0) {v += n;continue;}
796: 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);
797: for (j=0; j<n; j++) {
798: if (indexn[j] < 0) {v++; continue;}
799: 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);
800: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
801: }
802: }
803: return(0);
804: }
806: /* -----------------------------------------------------------------*/
810: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
811: {
812: Mat_SeqDense *a;
814: PetscInt *scols,i,j,nz,header[4];
815: int fd;
816: PetscMPIInt size;
817: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
818: PetscScalar *vals,*svals,*v,*w;
819: MPI_Comm comm;
822: PetscObjectGetComm((PetscObject)viewer,&comm);
823: MPI_Comm_size(comm,&size);
824: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
825: PetscViewerBinaryGetDescriptor(viewer,&fd);
826: PetscBinaryRead(fd,header,4,PETSC_INT);
827: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
828: M = header[1]; N = header[2]; nz = header[3];
830: /* set global size if not set already*/
831: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
832: MatSetSizes(newmat,M,N,M,N);
833: } else {
834: /* if sizes and type are already set, check if the vector global sizes are correct */
835: MatGetSize(newmat,&grows,&gcols);
836: 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);
837: }
838: a = (Mat_SeqDense*)newmat->data;
839: if (!a->user_alloc) {
840: MatSeqDenseSetPreallocation(newmat,NULL);
841: }
843: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
844: a = (Mat_SeqDense*)newmat->data;
845: v = a->v;
846: /* Allocate some temp space to read in the values and then flip them
847: from row major to column major */
848: PetscMalloc1((M*N > 0 ? M*N : 1),&w);
849: /* read in nonzero values */
850: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
851: /* now flip the values and store them in the matrix*/
852: for (j=0; j<N; j++) {
853: for (i=0; i<M; i++) {
854: *v++ =w[i*N+j];
855: }
856: }
857: PetscFree(w);
858: } else {
859: /* read row lengths */
860: PetscMalloc1((M+1),&rowlengths);
861: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
863: a = (Mat_SeqDense*)newmat->data;
864: v = a->v;
866: /* read column indices and nonzeros */
867: PetscMalloc1((nz+1),&scols);
868: cols = scols;
869: PetscBinaryRead(fd,cols,nz,PETSC_INT);
870: PetscMalloc1((nz+1),&svals);
871: vals = svals;
872: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
874: /* insert into matrix */
875: for (i=0; i<M; i++) {
876: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
877: svals += rowlengths[i]; scols += rowlengths[i];
878: }
879: PetscFree(vals);
880: PetscFree(cols);
881: PetscFree(rowlengths);
882: }
883: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
884: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
885: return(0);
886: }
890: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
891: {
892: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
893: PetscErrorCode ierr;
894: PetscInt i,j;
895: const char *name;
896: PetscScalar *v;
897: PetscViewerFormat format;
898: #if defined(PETSC_USE_COMPLEX)
899: PetscBool allreal = PETSC_TRUE;
900: #endif
903: PetscViewerGetFormat(viewer,&format);
904: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
905: return(0); /* do nothing for now */
906: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
907: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
908: for (i=0; i<A->rmap->n; i++) {
909: v = a->v + i;
910: PetscViewerASCIIPrintf(viewer,"row %D:",i);
911: for (j=0; j<A->cmap->n; j++) {
912: #if defined(PETSC_USE_COMPLEX)
913: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
914: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
915: } else if (PetscRealPart(*v)) {
916: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
917: }
918: #else
919: if (*v) {
920: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
921: }
922: #endif
923: v += a->lda;
924: }
925: PetscViewerASCIIPrintf(viewer,"\n");
926: }
927: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
928: } else {
929: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
930: #if defined(PETSC_USE_COMPLEX)
931: /* determine if matrix has all real values */
932: v = a->v;
933: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
934: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
935: }
936: #endif
937: if (format == PETSC_VIEWER_ASCII_MATLAB) {
938: PetscObjectGetName((PetscObject)A,&name);
939: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
940: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
941: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
942: }
944: for (i=0; i<A->rmap->n; i++) {
945: v = a->v + i;
946: for (j=0; j<A->cmap->n; j++) {
947: #if defined(PETSC_USE_COMPLEX)
948: if (allreal) {
949: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
950: } else {
951: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
952: }
953: #else
954: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
955: #endif
956: v += a->lda;
957: }
958: PetscViewerASCIIPrintf(viewer,"\n");
959: }
960: if (format == PETSC_VIEWER_ASCII_MATLAB) {
961: PetscViewerASCIIPrintf(viewer,"];\n");
962: }
963: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
964: }
965: PetscViewerFlush(viewer);
966: return(0);
967: }
971: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
972: {
973: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
974: PetscErrorCode ierr;
975: int fd;
976: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
977: PetscScalar *v,*anonz,*vals;
978: PetscViewerFormat format;
981: PetscViewerBinaryGetDescriptor(viewer,&fd);
983: PetscViewerGetFormat(viewer,&format);
984: if (format == PETSC_VIEWER_NATIVE) {
985: /* store the matrix as a dense matrix */
986: PetscMalloc1(4,&col_lens);
988: col_lens[0] = MAT_FILE_CLASSID;
989: col_lens[1] = m;
990: col_lens[2] = n;
991: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
993: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
994: PetscFree(col_lens);
996: /* write out matrix, by rows */
997: PetscMalloc1((m*n+1),&vals);
998: v = a->v;
999: for (j=0; j<n; j++) {
1000: for (i=0; i<m; i++) {
1001: vals[j + i*n] = *v++;
1002: }
1003: }
1004: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1005: PetscFree(vals);
1006: } else {
1007: PetscMalloc1((4+nz),&col_lens);
1009: col_lens[0] = MAT_FILE_CLASSID;
1010: col_lens[1] = m;
1011: col_lens[2] = n;
1012: col_lens[3] = nz;
1014: /* store lengths of each row and write (including header) to file */
1015: for (i=0; i<m; i++) col_lens[4+i] = n;
1016: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1018: /* Possibly should write in smaller increments, not whole matrix at once? */
1019: /* store column indices (zero start index) */
1020: ict = 0;
1021: for (i=0; i<m; i++) {
1022: for (j=0; j<n; j++) col_lens[ict++] = j;
1023: }
1024: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1025: PetscFree(col_lens);
1027: /* store nonzero values */
1028: PetscMalloc1((nz+1),&anonz);
1029: ict = 0;
1030: for (i=0; i<m; i++) {
1031: v = a->v + i;
1032: for (j=0; j<n; j++) {
1033: anonz[ict++] = *v; v += a->lda;
1034: }
1035: }
1036: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1037: PetscFree(anonz);
1038: }
1039: return(0);
1040: }
1042: #include <petscdraw.h>
1045: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1046: {
1047: Mat A = (Mat) Aa;
1048: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1049: PetscErrorCode ierr;
1050: PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j;
1051: PetscScalar *v = a->v;
1052: PetscViewer viewer;
1053: PetscDraw popup;
1054: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1055: PetscViewerFormat format;
1058: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1059: PetscViewerGetFormat(viewer,&format);
1060: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1062: /* Loop over matrix elements drawing boxes */
1063: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1064: /* Blue for negative and Red for positive */
1065: color = PETSC_DRAW_BLUE;
1066: for (j = 0; j < n; j++) {
1067: x_l = j;
1068: x_r = x_l + 1.0;
1069: for (i = 0; i < m; i++) {
1070: y_l = m - i - 1.0;
1071: y_r = y_l + 1.0;
1072: if (PetscRealPart(v[j*m+i]) > 0.) {
1073: color = PETSC_DRAW_RED;
1074: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1075: color = PETSC_DRAW_BLUE;
1076: } else {
1077: continue;
1078: }
1079: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1080: }
1081: }
1082: } else {
1083: /* use contour shading to indicate magnitude of values */
1084: /* first determine max of all nonzero values */
1085: for (i = 0; i < m*n; i++) {
1086: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1087: }
1088: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1089: PetscDrawGetPopup(draw,&popup);
1090: if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1091: for (j = 0; j < n; j++) {
1092: x_l = j;
1093: x_r = x_l + 1.0;
1094: for (i = 0; i < m; i++) {
1095: y_l = m - i - 1.0;
1096: y_r = y_l + 1.0;
1097: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1098: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1099: }
1100: }
1101: }
1102: return(0);
1103: }
1107: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1108: {
1109: PetscDraw draw;
1110: PetscBool isnull;
1111: PetscReal xr,yr,xl,yl,h,w;
1115: PetscViewerDrawGetDraw(viewer,0,&draw);
1116: PetscDrawIsNull(draw,&isnull);
1117: if (isnull) return(0);
1119: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1120: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1121: xr += w; yr += h; xl = -w; yl = -h;
1122: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1123: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1124: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1125: return(0);
1126: }
1130: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1131: {
1133: PetscBool iascii,isbinary,isdraw;
1136: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1137: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1138: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1140: if (iascii) {
1141: MatView_SeqDense_ASCII(A,viewer);
1142: } else if (isbinary) {
1143: MatView_SeqDense_Binary(A,viewer);
1144: } else if (isdraw) {
1145: MatView_SeqDense_Draw(A,viewer);
1146: }
1147: return(0);
1148: }
1152: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1153: {
1154: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1158: #if defined(PETSC_USE_LOG)
1159: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1160: #endif
1161: PetscFree(l->pivots);
1162: if (!l->user_alloc) {PetscFree(l->v);}
1163: PetscFree(mat->data);
1165: PetscObjectChangeTypeName((PetscObject)mat,0);
1166: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1167: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1168: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1169: #if defined(PETSC_HAVE_ELEMENTAL)
1170: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1171: #endif
1172: PetscObjectComposeFunction((PetscObject)mat,"MatGetFactor_petsc_C",NULL);
1173: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1174: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1175: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1176: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1177: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1178: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1179: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1180: return(0);
1181: }
1185: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1186: {
1187: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1189: PetscInt k,j,m,n,M;
1190: PetscScalar *v,tmp;
1193: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1194: if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1195: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1196: else {
1197: for (j=0; j<m; j++) {
1198: for (k=0; k<j; k++) {
1199: tmp = v[j + k*M];
1200: v[j + k*M] = v[k + j*M];
1201: v[k + j*M] = tmp;
1202: }
1203: }
1204: }
1205: } else { /* out-of-place transpose */
1206: Mat tmat;
1207: Mat_SeqDense *tmatd;
1208: PetscScalar *v2;
1210: if (reuse == MAT_INITIAL_MATRIX) {
1211: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1212: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1213: MatSetType(tmat,((PetscObject)A)->type_name);
1214: MatSeqDenseSetPreallocation(tmat,NULL);
1215: } else {
1216: tmat = *matout;
1217: }
1218: tmatd = (Mat_SeqDense*)tmat->data;
1219: v = mat->v; v2 = tmatd->v;
1220: for (j=0; j<n; j++) {
1221: for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
1222: }
1223: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1224: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1225: *matout = tmat;
1226: }
1227: return(0);
1228: }
1232: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1233: {
1234: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1235: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1236: PetscInt i,j;
1237: PetscScalar *v1,*v2;
1240: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1241: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1242: for (i=0; i<A1->rmap->n; i++) {
1243: v1 = mat1->v+i; v2 = mat2->v+i;
1244: for (j=0; j<A1->cmap->n; j++) {
1245: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1246: v1 += mat1->lda; v2 += mat2->lda;
1247: }
1248: }
1249: *flg = PETSC_TRUE;
1250: return(0);
1251: }
1255: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1256: {
1257: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1259: PetscInt i,n,len;
1260: PetscScalar *x,zero = 0.0;
1263: VecSet(v,zero);
1264: VecGetSize(v,&n);
1265: VecGetArray(v,&x);
1266: len = PetscMin(A->rmap->n,A->cmap->n);
1267: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1268: for (i=0; i<len; i++) {
1269: x[i] = mat->v[i*mat->lda + i];
1270: }
1271: VecRestoreArray(v,&x);
1272: return(0);
1273: }
1277: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1278: {
1279: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1280: PetscScalar *l,*r,x,*v;
1282: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1285: if (ll) {
1286: VecGetSize(ll,&m);
1287: VecGetArray(ll,&l);
1288: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1289: for (i=0; i<m; i++) {
1290: x = l[i];
1291: v = mat->v + i;
1292: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1293: }
1294: VecRestoreArray(ll,&l);
1295: PetscLogFlops(n*m);
1296: }
1297: if (rr) {
1298: VecGetSize(rr,&n);
1299: VecGetArray(rr,&r);
1300: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1301: for (i=0; i<n; i++) {
1302: x = r[i];
1303: v = mat->v + i*m;
1304: for (j=0; j<m; j++) (*v++) *= x;
1305: }
1306: VecRestoreArray(rr,&r);
1307: PetscLogFlops(n*m);
1308: }
1309: return(0);
1310: }
1314: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1315: {
1316: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1317: PetscScalar *v = mat->v;
1318: PetscReal sum = 0.0;
1319: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1323: if (type == NORM_FROBENIUS) {
1324: if (lda>m) {
1325: for (j=0; j<A->cmap->n; j++) {
1326: v = mat->v+j*lda;
1327: for (i=0; i<m; i++) {
1328: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1329: }
1330: }
1331: } else {
1332: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1333: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1334: }
1335: }
1336: *nrm = PetscSqrtReal(sum);
1337: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1338: } else if (type == NORM_1) {
1339: *nrm = 0.0;
1340: for (j=0; j<A->cmap->n; j++) {
1341: v = mat->v + j*mat->lda;
1342: sum = 0.0;
1343: for (i=0; i<A->rmap->n; i++) {
1344: sum += PetscAbsScalar(*v); v++;
1345: }
1346: if (sum > *nrm) *nrm = sum;
1347: }
1348: PetscLogFlops(A->cmap->n*A->rmap->n);
1349: } else if (type == NORM_INFINITY) {
1350: *nrm = 0.0;
1351: for (j=0; j<A->rmap->n; j++) {
1352: v = mat->v + j;
1353: sum = 0.0;
1354: for (i=0; i<A->cmap->n; i++) {
1355: sum += PetscAbsScalar(*v); v += mat->lda;
1356: }
1357: if (sum > *nrm) *nrm = sum;
1358: }
1359: PetscLogFlops(A->cmap->n*A->rmap->n);
1360: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1361: return(0);
1362: }
1366: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1367: {
1368: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1372: switch (op) {
1373: case MAT_ROW_ORIENTED:
1374: aij->roworiented = flg;
1375: break;
1376: case MAT_NEW_NONZERO_LOCATIONS:
1377: case MAT_NEW_NONZERO_LOCATION_ERR:
1378: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1379: case MAT_NEW_DIAGONALS:
1380: case MAT_KEEP_NONZERO_PATTERN:
1381: case MAT_IGNORE_OFF_PROC_ENTRIES:
1382: case MAT_USE_HASH_TABLE:
1383: case MAT_IGNORE_LOWER_TRIANGULAR:
1384: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1385: break;
1386: case MAT_SPD:
1387: case MAT_SYMMETRIC:
1388: case MAT_STRUCTURALLY_SYMMETRIC:
1389: case MAT_HERMITIAN:
1390: case MAT_SYMMETRY_ETERNAL:
1391: /* These options are handled directly by MatSetOption() */
1392: break;
1393: default:
1394: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1395: }
1396: return(0);
1397: }
1401: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1402: {
1403: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1405: PetscInt lda=l->lda,m=A->rmap->n,j;
1408: if (lda>m) {
1409: for (j=0; j<A->cmap->n; j++) {
1410: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1411: }
1412: } else {
1413: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1414: }
1415: return(0);
1416: }
1420: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1421: {
1422: PetscErrorCode ierr;
1423: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1424: PetscInt m = l->lda, n = A->cmap->n, i,j;
1425: PetscScalar *slot,*bb;
1426: const PetscScalar *xx;
1429: #if defined(PETSC_USE_DEBUG)
1430: for (i=0; i<N; i++) {
1431: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1432: 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);
1433: }
1434: #endif
1436: /* fix right hand side if needed */
1437: if (x && b) {
1438: VecGetArrayRead(x,&xx);
1439: VecGetArray(b,&bb);
1440: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1441: VecRestoreArrayRead(x,&xx);
1442: VecRestoreArray(b,&bb);
1443: }
1445: for (i=0; i<N; i++) {
1446: slot = l->v + rows[i];
1447: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1448: }
1449: if (diag != 0.0) {
1450: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1451: for (i=0; i<N; i++) {
1452: slot = l->v + (m+1)*rows[i];
1453: *slot = diag;
1454: }
1455: }
1456: return(0);
1457: }
1461: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1462: {
1463: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1466: 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");
1467: *array = mat->v;
1468: return(0);
1469: }
1473: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1474: {
1476: *array = 0; /* user cannot accidently use the array later */
1477: return(0);
1478: }
1482: /*@C
1483: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1485: Not Collective
1487: Input Parameter:
1488: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1490: Output Parameter:
1491: . array - pointer to the data
1493: Level: intermediate
1495: .seealso: MatDenseRestoreArray()
1496: @*/
1497: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1498: {
1502: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1503: return(0);
1504: }
1508: /*@C
1509: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1511: Not Collective
1513: Input Parameters:
1514: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1515: . array - pointer to the data
1517: Level: intermediate
1519: .seealso: MatDenseGetArray()
1520: @*/
1521: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1522: {
1526: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1527: return(0);
1528: }
1532: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1533: {
1534: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1536: PetscInt i,j,nrows,ncols;
1537: const PetscInt *irow,*icol;
1538: PetscScalar *av,*bv,*v = mat->v;
1539: Mat newmat;
1542: ISGetIndices(isrow,&irow);
1543: ISGetIndices(iscol,&icol);
1544: ISGetLocalSize(isrow,&nrows);
1545: ISGetLocalSize(iscol,&ncols);
1547: /* Check submatrixcall */
1548: if (scall == MAT_REUSE_MATRIX) {
1549: PetscInt n_cols,n_rows;
1550: MatGetSize(*B,&n_rows,&n_cols);
1551: if (n_rows != nrows || n_cols != ncols) {
1552: /* resize the result matrix to match number of requested rows/columns */
1553: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1554: }
1555: newmat = *B;
1556: } else {
1557: /* Create and fill new matrix */
1558: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1559: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1560: MatSetType(newmat,((PetscObject)A)->type_name);
1561: MatSeqDenseSetPreallocation(newmat,NULL);
1562: }
1564: /* Now extract the data pointers and do the copy,column at a time */
1565: bv = ((Mat_SeqDense*)newmat->data)->v;
1567: for (i=0; i<ncols; i++) {
1568: av = v + mat->lda*icol[i];
1569: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1570: }
1572: /* Assemble the matrices so that the correct flags are set */
1573: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1574: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1576: /* Free work space */
1577: ISRestoreIndices(isrow,&irow);
1578: ISRestoreIndices(iscol,&icol);
1579: *B = newmat;
1580: return(0);
1581: }
1585: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1586: {
1588: PetscInt i;
1591: if (scall == MAT_INITIAL_MATRIX) {
1592: PetscMalloc1((n+1),B);
1593: }
1595: for (i=0; i<n; i++) {
1596: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1597: }
1598: return(0);
1599: }
1603: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1604: {
1606: return(0);
1607: }
1611: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1612: {
1614: return(0);
1615: }
1619: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1620: {
1621: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1623: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1626: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1627: if (A->ops->copy != B->ops->copy) {
1628: MatCopy_Basic(A,B,str);
1629: return(0);
1630: }
1631: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1632: if (lda1>m || lda2>m) {
1633: for (j=0; j<n; j++) {
1634: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1635: }
1636: } else {
1637: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1638: }
1639: return(0);
1640: }
1644: PetscErrorCode MatSetUp_SeqDense(Mat A)
1645: {
1649: MatSeqDenseSetPreallocation(A,0);
1650: return(0);
1651: }
1655: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1656: {
1657: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1658: PetscInt i,nz = A->rmap->n*A->cmap->n;
1659: PetscScalar *aa = a->v;
1662: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1663: return(0);
1664: }
1668: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1669: {
1670: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1671: PetscInt i,nz = A->rmap->n*A->cmap->n;
1672: PetscScalar *aa = a->v;
1675: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1676: return(0);
1677: }
1681: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1682: {
1683: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1684: PetscInt i,nz = A->rmap->n*A->cmap->n;
1685: PetscScalar *aa = a->v;
1688: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1689: return(0);
1690: }
1692: /* ----------------------------------------------------------------*/
1695: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1696: {
1700: if (scall == MAT_INITIAL_MATRIX) {
1701: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1702: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1703: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1704: }
1705: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1706: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1707: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1708: return(0);
1709: }
1713: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1714: {
1716: PetscInt m=A->rmap->n,n=B->cmap->n;
1717: Mat Cmat;
1720: 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);
1721: MatCreate(PETSC_COMM_SELF,&Cmat);
1722: MatSetSizes(Cmat,m,n,m,n);
1723: MatSetType(Cmat,MATSEQDENSE);
1724: MatSeqDenseSetPreallocation(Cmat,NULL);
1726: *C = Cmat;
1727: return(0);
1728: }
1732: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1733: {
1734: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1735: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1736: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1737: PetscBLASInt m,n,k;
1738: PetscScalar _DOne=1.0,_DZero=0.0;
1742: PetscBLASIntCast(A->rmap->n,&m);
1743: PetscBLASIntCast(B->cmap->n,&n);
1744: PetscBLASIntCast(A->cmap->n,&k);
1745: PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1746: return(0);
1747: }
1751: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1752: {
1756: if (scall == MAT_INITIAL_MATRIX) {
1757: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1758: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1759: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1760: }
1761: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1762: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1763: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1764: return(0);
1765: }
1769: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1770: {
1772: PetscInt m=A->cmap->n,n=B->cmap->n;
1773: Mat Cmat;
1776: 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);
1777: MatCreate(PETSC_COMM_SELF,&Cmat);
1778: MatSetSizes(Cmat,m,n,m,n);
1779: MatSetType(Cmat,MATSEQDENSE);
1780: MatSeqDenseSetPreallocation(Cmat,NULL);
1782: Cmat->assembled = PETSC_TRUE;
1784: *C = Cmat;
1785: return(0);
1786: }
1790: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1791: {
1792: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1793: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1794: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1795: PetscBLASInt m,n,k;
1796: PetscScalar _DOne=1.0,_DZero=0.0;
1800: PetscBLASIntCast(A->cmap->n,&m);
1801: PetscBLASIntCast(B->cmap->n,&n);
1802: PetscBLASIntCast(A->rmap->n,&k);
1803: /*
1804: Note the m and n arguments below are the number rows and columns of A', not A!
1805: */
1806: PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1807: return(0);
1808: }
1812: PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1813: {
1814: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1816: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1817: PetscScalar *x;
1818: MatScalar *aa = a->v;
1821: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1823: VecSet(v,0.0);
1824: VecGetArray(v,&x);
1825: VecGetLocalSize(v,&p);
1826: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1827: for (i=0; i<m; i++) {
1828: x[i] = aa[i]; if (idx) idx[i] = 0;
1829: for (j=1; j<n; j++) {
1830: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1831: }
1832: }
1833: VecRestoreArray(v,&x);
1834: return(0);
1835: }
1839: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1840: {
1841: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1843: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1844: PetscScalar *x;
1845: PetscReal atmp;
1846: MatScalar *aa = a->v;
1849: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1851: VecSet(v,0.0);
1852: VecGetArray(v,&x);
1853: VecGetLocalSize(v,&p);
1854: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1855: for (i=0; i<m; i++) {
1856: x[i] = PetscAbsScalar(aa[i]);
1857: for (j=1; j<n; j++) {
1858: atmp = PetscAbsScalar(aa[i+m*j]);
1859: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1860: }
1861: }
1862: VecRestoreArray(v,&x);
1863: return(0);
1864: }
1868: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1869: {
1870: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1872: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1873: PetscScalar *x;
1874: MatScalar *aa = a->v;
1877: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1879: VecSet(v,0.0);
1880: VecGetArray(v,&x);
1881: VecGetLocalSize(v,&p);
1882: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1883: for (i=0; i<m; i++) {
1884: x[i] = aa[i]; if (idx) idx[i] = 0;
1885: for (j=1; j<n; j++) {
1886: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1887: }
1888: }
1889: VecRestoreArray(v,&x);
1890: return(0);
1891: }
1895: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1896: {
1897: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1899: PetscScalar *x;
1902: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1904: VecGetArray(v,&x);
1905: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1906: VecRestoreArray(v,&x);
1907: return(0);
1908: }
1913: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1914: {
1916: PetscInt i,j,m,n;
1917: PetscScalar *a;
1920: MatGetSize(A,&m,&n);
1921: PetscMemzero(norms,n*sizeof(PetscReal));
1922: MatDenseGetArray(A,&a);
1923: if (type == NORM_2) {
1924: for (i=0; i<n; i++) {
1925: for (j=0; j<m; j++) {
1926: norms[i] += PetscAbsScalar(a[j]*a[j]);
1927: }
1928: a += m;
1929: }
1930: } else if (type == NORM_1) {
1931: for (i=0; i<n; i++) {
1932: for (j=0; j<m; j++) {
1933: norms[i] += PetscAbsScalar(a[j]);
1934: }
1935: a += m;
1936: }
1937: } else if (type == NORM_INFINITY) {
1938: for (i=0; i<n; i++) {
1939: for (j=0; j<m; j++) {
1940: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1941: }
1942: a += m;
1943: }
1944: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1945: MatDenseRestoreArray(A,&a);
1946: if (type == NORM_2) {
1947: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1948: }
1949: return(0);
1950: }
1954: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
1955: {
1957: PetscScalar *a;
1958: PetscInt m,n,i;
1961: MatGetSize(x,&m,&n);
1962: MatDenseGetArray(x,&a);
1963: for (i=0; i<m*n; i++) {
1964: PetscRandomGetValue(rctx,a+i);
1965: }
1966: MatDenseRestoreArray(x,&a);
1967: return(0);
1968: }
1971: /* -------------------------------------------------------------------*/
1972: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
1973: MatGetRow_SeqDense,
1974: MatRestoreRow_SeqDense,
1975: MatMult_SeqDense,
1976: /* 4*/ MatMultAdd_SeqDense,
1977: MatMultTranspose_SeqDense,
1978: MatMultTransposeAdd_SeqDense,
1979: 0,
1980: 0,
1981: 0,
1982: /* 10*/ 0,
1983: MatLUFactor_SeqDense,
1984: MatCholeskyFactor_SeqDense,
1985: MatSOR_SeqDense,
1986: MatTranspose_SeqDense,
1987: /* 15*/ MatGetInfo_SeqDense,
1988: MatEqual_SeqDense,
1989: MatGetDiagonal_SeqDense,
1990: MatDiagonalScale_SeqDense,
1991: MatNorm_SeqDense,
1992: /* 20*/ MatAssemblyBegin_SeqDense,
1993: MatAssemblyEnd_SeqDense,
1994: MatSetOption_SeqDense,
1995: MatZeroEntries_SeqDense,
1996: /* 24*/ MatZeroRows_SeqDense,
1997: 0,
1998: 0,
1999: 0,
2000: 0,
2001: /* 29*/ MatSetUp_SeqDense,
2002: 0,
2003: 0,
2004: 0,
2005: 0,
2006: /* 34*/ MatDuplicate_SeqDense,
2007: 0,
2008: 0,
2009: 0,
2010: 0,
2011: /* 39*/ MatAXPY_SeqDense,
2012: MatGetSubMatrices_SeqDense,
2013: 0,
2014: MatGetValues_SeqDense,
2015: MatCopy_SeqDense,
2016: /* 44*/ MatGetRowMax_SeqDense,
2017: MatScale_SeqDense,
2018: 0,
2019: 0,
2020: 0,
2021: /* 49*/ MatSetRandom_SeqDense,
2022: 0,
2023: 0,
2024: 0,
2025: 0,
2026: /* 54*/ 0,
2027: 0,
2028: 0,
2029: 0,
2030: 0,
2031: /* 59*/ 0,
2032: MatDestroy_SeqDense,
2033: MatView_SeqDense,
2034: 0,
2035: 0,
2036: /* 64*/ 0,
2037: 0,
2038: 0,
2039: 0,
2040: 0,
2041: /* 69*/ MatGetRowMaxAbs_SeqDense,
2042: 0,
2043: 0,
2044: 0,
2045: 0,
2046: /* 74*/ 0,
2047: 0,
2048: 0,
2049: 0,
2050: 0,
2051: /* 79*/ 0,
2052: 0,
2053: 0,
2054: 0,
2055: /* 83*/ MatLoad_SeqDense,
2056: 0,
2057: MatIsHermitian_SeqDense,
2058: 0,
2059: 0,
2060: 0,
2061: /* 89*/ MatMatMult_SeqDense_SeqDense,
2062: MatMatMultSymbolic_SeqDense_SeqDense,
2063: MatMatMultNumeric_SeqDense_SeqDense,
2064: 0,
2065: 0,
2066: /* 94*/ 0,
2067: 0,
2068: 0,
2069: 0,
2070: 0,
2071: /* 99*/ 0,
2072: 0,
2073: 0,
2074: MatConjugate_SeqDense,
2075: 0,
2076: /*104*/ 0,
2077: MatRealPart_SeqDense,
2078: MatImaginaryPart_SeqDense,
2079: 0,
2080: 0,
2081: /*109*/ MatMatSolve_SeqDense,
2082: 0,
2083: MatGetRowMin_SeqDense,
2084: MatGetColumnVector_SeqDense,
2085: 0,
2086: /*114*/ 0,
2087: 0,
2088: 0,
2089: 0,
2090: 0,
2091: /*119*/ 0,
2092: 0,
2093: 0,
2094: 0,
2095: 0,
2096: /*124*/ 0,
2097: MatGetColumnNorms_SeqDense,
2098: 0,
2099: 0,
2100: 0,
2101: /*129*/ 0,
2102: MatTransposeMatMult_SeqDense_SeqDense,
2103: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2104: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2105: 0,
2106: /*134*/ 0,
2107: 0,
2108: 0,
2109: 0,
2110: 0,
2111: /*139*/ 0,
2112: 0,
2113: 0
2114: };
2118: /*@C
2119: MatCreateSeqDense - Creates a sequential dense matrix that
2120: is stored in column major order (the usual Fortran 77 manner). Many
2121: of the matrix operations use the BLAS and LAPACK routines.
2123: Collective on MPI_Comm
2125: Input Parameters:
2126: + comm - MPI communicator, set to PETSC_COMM_SELF
2127: . m - number of rows
2128: . n - number of columns
2129: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2130: to control all matrix memory allocation.
2132: Output Parameter:
2133: . A - the matrix
2135: Notes:
2136: The data input variable is intended primarily for Fortran programmers
2137: who wish to allocate their own matrix memory space. Most users should
2138: set data=NULL.
2140: Level: intermediate
2142: .keywords: dense, matrix, LAPACK, BLAS
2144: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2145: @*/
2146: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2147: {
2151: MatCreate(comm,A);
2152: MatSetSizes(*A,m,n,m,n);
2153: MatSetType(*A,MATSEQDENSE);
2154: MatSeqDenseSetPreallocation(*A,data);
2155: return(0);
2156: }
2160: /*@C
2161: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2163: Collective on MPI_Comm
2165: Input Parameters:
2166: + B - the matrix
2167: - data - the array (or NULL)
2169: Notes:
2170: The data input variable is intended primarily for Fortran programmers
2171: who wish to allocate their own matrix memory space. Most users should
2172: need not call this routine.
2174: Level: intermediate
2176: .keywords: dense, matrix, LAPACK, BLAS
2178: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2180: @*/
2181: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2182: {
2186: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2187: return(0);
2188: }
2192: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2193: {
2194: Mat_SeqDense *b;
2198: B->preallocated = PETSC_TRUE;
2200: PetscLayoutSetUp(B->rmap);
2201: PetscLayoutSetUp(B->cmap);
2203: b = (Mat_SeqDense*)B->data;
2204: b->Mmax = B->rmap->n;
2205: b->Nmax = B->cmap->n;
2206: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2208: if (!data) { /* petsc-allocated storage */
2209: if (!b->user_alloc) { PetscFree(b->v); }
2210: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2211: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2213: b->user_alloc = PETSC_FALSE;
2214: } else { /* user-allocated storage */
2215: if (!b->user_alloc) { PetscFree(b->v); }
2216: b->v = data;
2217: b->user_alloc = PETSC_TRUE;
2218: }
2219: B->assembled = PETSC_TRUE;
2220: return(0);
2221: }
2225: PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2226: {
2228: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for requested operation");
2229: return(0);
2230: }
2234: /*@C
2235: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2237: Input parameter:
2238: + A - the matrix
2239: - lda - the leading dimension
2241: Notes:
2242: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2243: it asserts that the preallocation has a leading dimension (the LDA parameter
2244: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2246: Level: intermediate
2248: .keywords: dense, matrix, LAPACK, BLAS
2250: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2252: @*/
2253: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2254: {
2255: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2258: 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);
2259: b->lda = lda;
2260: b->changelda = PETSC_FALSE;
2261: b->Mmax = PetscMax(b->Mmax,lda);
2262: return(0);
2263: }
2265: /*MC
2266: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2268: Options Database Keys:
2269: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2271: Level: beginner
2273: .seealso: MatCreateSeqDense()
2275: M*/
2279: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2280: {
2281: Mat_SeqDense *b;
2283: PetscMPIInt size;
2286: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2287: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2289: PetscNewLog(B,&b);
2290: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2291: B->data = (void*)b;
2293: b->pivots = 0;
2294: b->roworiented = PETSC_TRUE;
2295: b->v = 0;
2296: b->changelda = PETSC_FALSE;
2298: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2299: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2300: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2301: #if defined(PETSC_HAVE_ELEMENTAL)
2302: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2303: #endif
2304: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);
2305: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2306: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2307: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2308: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2309: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2310: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2311: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2312: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2313: return(0);
2314: }