Actual source code: dense.c
petsc-3.4.5 2014-06-29
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;
19: PetscInt *rows;
20: MatScalar *aa = a->v;
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: MatSeqAIJSetPreallocation(B,A->cmap->n,NULL);
28: PetscMalloc(A->rmap->n*sizeof(PetscInt),&rows);
29: for (i=0; i<A->rmap->n; i++) rows[i] = i;
31: for (i=0; i<A->cmap->n; i++) {
32: MatSetValues(B,A->rmap->n,rows,1,&i,aa,INSERT_VALUES);
33: aa += a->lda;
34: }
35: PetscFree(rows);
36: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
37: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
39: if (reuse == MAT_REUSE_MATRIX) {
40: MatHeaderReplace(A,B);
41: } else {
42: *newmat = B;
43: }
44: return(0);
45: }
49: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
50: {
51: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
52: PetscScalar oalpha = alpha;
53: PetscInt j;
54: PetscBLASInt N,m,ldax,lday,one = 1;
58: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
59: PetscBLASIntCast(X->rmap->n,&m);
60: PetscBLASIntCast(x->lda,&ldax);
61: PetscBLASIntCast(y->lda,&lday);
62: if (ldax>m || lday>m) {
63: for (j=0; j<X->cmap->n; j++) {
64: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
65: }
66: } else {
67: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
68: }
69: PetscLogFlops(PetscMax(2*N-1,0));
70: return(0);
71: }
75: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
76: {
77: PetscInt N = A->rmap->n*A->cmap->n;
80: info->block_size = 1.0;
81: info->nz_allocated = (double)N;
82: info->nz_used = (double)N;
83: info->nz_unneeded = (double)0;
84: info->assemblies = (double)A->num_ass;
85: info->mallocs = 0;
86: info->memory = ((PetscObject)A)->mem;
87: info->fill_ratio_given = 0;
88: info->fill_ratio_needed = 0;
89: info->factor_mallocs = 0;
90: return(0);
91: }
95: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
96: {
97: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
98: PetscScalar oalpha = alpha;
100: PetscBLASInt one = 1,j,nz,lda;
103: PetscBLASIntCast(a->lda,&lda);
104: if (lda>A->rmap->n) {
105: PetscBLASIntCast(A->rmap->n,&nz);
106: for (j=0; j<A->cmap->n; j++) {
107: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
108: }
109: } else {
110: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
111: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
112: }
113: PetscLogFlops(nz);
114: return(0);
115: }
119: PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
120: {
121: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
122: PetscInt i,j,m = A->rmap->n,N;
123: PetscScalar *v = a->v;
126: *fl = PETSC_FALSE;
127: if (A->rmap->n != A->cmap->n) return(0);
128: N = a->lda;
130: for (i=0; i<m; i++) {
131: for (j=i+1; j<m; j++) {
132: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
133: }
134: }
135: *fl = PETSC_TRUE;
136: return(0);
137: }
141: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
142: {
143: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
145: PetscInt lda = (PetscInt)mat->lda,j,m;
148: PetscLayoutReference(A->rmap,&newi->rmap);
149: PetscLayoutReference(A->cmap,&newi->cmap);
150: MatSeqDenseSetPreallocation(newi,NULL);
151: if (cpvalues == MAT_COPY_VALUES) {
152: l = (Mat_SeqDense*)newi->data;
153: if (lda>A->rmap->n) {
154: m = A->rmap->n;
155: for (j=0; j<A->cmap->n; j++) {
156: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
157: }
158: } else {
159: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
160: }
161: }
162: newi->assembled = PETSC_TRUE;
163: return(0);
164: }
168: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
169: {
173: MatCreate(PetscObjectComm((PetscObject)A),newmat);
174: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
175: MatSetType(*newmat,((PetscObject)A)->type_name);
176: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
177: return(0);
178: }
181: extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
185: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
186: {
187: MatFactorInfo info;
191: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
192: MatLUFactor_SeqDense(fact,0,0,&info);
193: return(0);
194: }
198: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
199: {
200: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
202: PetscScalar *x,*y;
203: PetscBLASInt one = 1,info,m;
206: PetscBLASIntCast(A->rmap->n,&m);
207: VecGetArray(xx,&x);
208: VecGetArray(yy,&y);
209: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
210: if (A->factortype == MAT_FACTOR_LU) {
211: #if defined(PETSC_MISSING_LAPACK_GETRS)
212: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
213: #else
214: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
215: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
216: #endif
217: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
218: #if defined(PETSC_MISSING_LAPACK_POTRS)
219: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
220: #else
221: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
222: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
223: #endif
224: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
225: VecRestoreArray(xx,&x);
226: VecRestoreArray(yy,&y);
227: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
228: return(0);
229: }
233: PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
234: {
235: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
237: PetscScalar *b,*x;
238: PetscInt n;
239: PetscBLASInt nrhs,info,m;
240: PetscBool flg;
243: PetscBLASIntCast(A->rmap->n,&m);
244: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
245: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
246: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
247: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
249: MatGetSize(B,NULL,&n);
250: PetscBLASIntCast(n,&nrhs);
251: MatDenseGetArray(B,&b);
252: MatDenseGetArray(X,&x);
254: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
256: if (A->factortype == MAT_FACTOR_LU) {
257: #if defined(PETSC_MISSING_LAPACK_GETRS)
258: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
259: #else
260: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
261: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
262: #endif
263: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
264: #if defined(PETSC_MISSING_LAPACK_POTRS)
265: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
266: #else
267: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
268: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
269: #endif
270: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
272: MatDenseRestoreArray(B,&b);
273: MatDenseRestoreArray(X,&x);
274: PetscLogFlops(nrhs*(2.0*m*m - m));
275: return(0);
276: }
280: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
281: {
282: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
284: PetscScalar *x,*y;
285: PetscBLASInt one = 1,info,m;
288: PetscBLASIntCast(A->rmap->n,&m);
289: VecGetArray(xx,&x);
290: VecGetArray(yy,&y);
291: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
292: /* assume if pivots exist then use LU; else Cholesky */
293: if (mat->pivots) {
294: #if defined(PETSC_MISSING_LAPACK_GETRS)
295: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
296: #else
297: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
298: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
299: #endif
300: } else {
301: #if defined(PETSC_MISSING_LAPACK_POTRS)
302: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
303: #else
304: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
305: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
306: #endif
307: }
308: VecRestoreArray(xx,&x);
309: VecRestoreArray(yy,&y);
310: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
311: return(0);
312: }
316: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
317: {
318: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
320: PetscScalar *x,*y,sone = 1.0;
321: Vec tmp = 0;
322: PetscBLASInt one = 1,info,m;
325: PetscBLASIntCast(A->rmap->n,&m);
326: VecGetArray(xx,&x);
327: VecGetArray(yy,&y);
328: if (!A->rmap->n || !A->cmap->n) return(0);
329: if (yy == zz) {
330: VecDuplicate(yy,&tmp);
331: PetscLogObjectParent(A,tmp);
332: VecCopy(yy,tmp);
333: }
334: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
335: /* assume if pivots exist then use LU; else Cholesky */
336: if (mat->pivots) {
337: #if defined(PETSC_MISSING_LAPACK_GETRS)
338: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
339: #else
340: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
341: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
342: #endif
343: } else {
344: #if defined(PETSC_MISSING_LAPACK_POTRS)
345: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
346: #else
347: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
348: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
349: #endif
350: }
351: if (tmp) {
352: VecAXPY(yy,sone,tmp);
353: VecDestroy(&tmp);
354: } else {
355: VecAXPY(yy,sone,zz);
356: }
357: VecRestoreArray(xx,&x);
358: VecRestoreArray(yy,&y);
359: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
360: return(0);
361: }
365: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
366: {
367: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
369: PetscScalar *x,*y,sone = 1.0;
370: Vec tmp;
371: PetscBLASInt one = 1,info,m;
374: PetscBLASIntCast(A->rmap->n,&m);
375: if (!A->rmap->n || !A->cmap->n) return(0);
376: VecGetArray(xx,&x);
377: VecGetArray(yy,&y);
378: if (yy == zz) {
379: VecDuplicate(yy,&tmp);
380: PetscLogObjectParent(A,tmp);
381: VecCopy(yy,tmp);
382: }
383: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
384: /* assume if pivots exist then use LU; else Cholesky */
385: if (mat->pivots) {
386: #if defined(PETSC_MISSING_LAPACK_GETRS)
387: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
388: #else
389: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
390: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
391: #endif
392: } else {
393: #if defined(PETSC_MISSING_LAPACK_POTRS)
394: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
395: #else
396: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
397: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
398: #endif
399: }
400: if (tmp) {
401: VecAXPY(yy,sone,tmp);
402: VecDestroy(&tmp);
403: } else {
404: VecAXPY(yy,sone,zz);
405: }
406: VecRestoreArray(xx,&x);
407: VecRestoreArray(yy,&y);
408: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
409: return(0);
410: }
412: /* ---------------------------------------------------------------*/
413: /* COMMENT: I have chosen to hide row permutation in the pivots,
414: rather than put it in the Mat->row slot.*/
417: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
418: {
419: #if defined(PETSC_MISSING_LAPACK_GETRF)
421: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
422: #else
423: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
425: PetscBLASInt n,m,info;
428: PetscBLASIntCast(A->cmap->n,&n);
429: PetscBLASIntCast(A->rmap->n,&m);
430: if (!mat->pivots) {
431: PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);
432: PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));
433: }
434: if (!A->rmap->n || !A->cmap->n) return(0);
435: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
436: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
437: PetscFPTrapPop();
439: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
440: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
441: A->ops->solve = MatSolve_SeqDense;
442: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
443: A->ops->solveadd = MatSolveAdd_SeqDense;
444: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
445: A->factortype = MAT_FACTOR_LU;
447: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
448: #endif
449: return(0);
450: }
454: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
455: {
456: #if defined(PETSC_MISSING_LAPACK_POTRF)
458: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
459: #else
460: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
462: PetscBLASInt info,n;
465: PetscBLASIntCast(A->cmap->n,&n);
466: PetscFree(mat->pivots);
468: if (!A->rmap->n || !A->cmap->n) return(0);
469: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
470: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
471: A->ops->solve = MatSolve_SeqDense;
472: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
473: A->ops->solveadd = MatSolveAdd_SeqDense;
474: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
475: A->factortype = MAT_FACTOR_CHOLESKY;
477: PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
478: #endif
479: return(0);
480: }
485: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
486: {
488: MatFactorInfo info;
491: info.fill = 1.0;
493: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
494: MatCholeskyFactor_SeqDense(fact,0,&info);
495: return(0);
496: }
500: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
501: {
503: fact->assembled = PETSC_TRUE;
504: fact->preallocated = PETSC_TRUE;
505: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
506: return(0);
507: }
511: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
512: {
514: fact->preallocated = PETSC_TRUE;
515: fact->assembled = PETSC_TRUE;
516: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
517: return(0);
518: }
522: PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
523: {
527: MatCreate(PetscObjectComm((PetscObject)A),fact);
528: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
529: MatSetType(*fact,((PetscObject)A)->type_name);
530: if (ftype == MAT_FACTOR_LU) {
531: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
532: } else {
533: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
534: }
535: (*fact)->factortype = ftype;
536: return(0);
537: }
539: /* ------------------------------------------------------------------*/
542: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
543: {
544: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
545: PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt;
547: PetscInt m = A->rmap->n,i;
548: PetscBLASInt o = 1,bm;
551: PetscBLASIntCast(m,&bm);
552: if (flag & SOR_ZERO_INITIAL_GUESS) {
553: /* this is a hack fix, should have another version without the second BLASdot */
554: VecSet(xx,zero);
555: }
556: VecGetArray(xx,&x);
557: VecGetArray(bb,&b);
558: its = its*lits;
559: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
560: while (its--) {
561: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
562: for (i=0; i<m; i++) {
563: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
564: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
565: }
566: }
567: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
568: for (i=m-1; i>=0; i--) {
569: PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o));
570: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
571: }
572: }
573: }
574: VecRestoreArray(bb,&b);
575: VecRestoreArray(xx,&x);
576: return(0);
577: }
579: /* -----------------------------------------------------------------*/
582: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
583: {
584: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
585: PetscScalar *v = mat->v,*x,*y;
587: PetscBLASInt m, n,_One=1;
588: PetscScalar _DOne=1.0,_DZero=0.0;
591: PetscBLASIntCast(A->rmap->n,&m);
592: PetscBLASIntCast(A->cmap->n,&n);
593: if (!A->rmap->n || !A->cmap->n) return(0);
594: VecGetArray(xx,&x);
595: VecGetArray(yy,&y);
596: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
597: VecRestoreArray(xx,&x);
598: VecRestoreArray(yy,&y);
599: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
600: return(0);
601: }
605: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
606: {
607: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
608: PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
610: PetscBLASInt m, n, _One=1;
613: PetscBLASIntCast(A->rmap->n,&m);
614: PetscBLASIntCast(A->cmap->n,&n);
615: if (!A->rmap->n || !A->cmap->n) return(0);
616: VecGetArray(xx,&x);
617: VecGetArray(yy,&y);
618: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
619: VecRestoreArray(xx,&x);
620: VecRestoreArray(yy,&y);
621: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
622: return(0);
623: }
627: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
628: {
629: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
630: PetscScalar *v = mat->v,*x,*y,_DOne=1.0;
632: PetscBLASInt m, n, _One=1;
635: PetscBLASIntCast(A->rmap->n,&m);
636: PetscBLASIntCast(A->cmap->n,&n);
637: if (!A->rmap->n || !A->cmap->n) return(0);
638: if (zz != yy) {VecCopy(zz,yy);}
639: VecGetArray(xx,&x);
640: VecGetArray(yy,&y);
641: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
642: VecRestoreArray(xx,&x);
643: VecRestoreArray(yy,&y);
644: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
645: return(0);
646: }
650: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
651: {
652: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
653: PetscScalar *v = mat->v,*x,*y;
655: PetscBLASInt m, n, _One=1;
656: PetscScalar _DOne=1.0;
659: PetscBLASIntCast(A->rmap->n,&m);
660: PetscBLASIntCast(A->cmap->n,&n);
661: if (!A->rmap->n || !A->cmap->n) return(0);
662: if (zz != yy) {VecCopy(zz,yy);}
663: VecGetArray(xx,&x);
664: VecGetArray(yy,&y);
665: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
666: VecRestoreArray(xx,&x);
667: VecRestoreArray(yy,&y);
668: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
669: return(0);
670: }
672: /* -----------------------------------------------------------------*/
675: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
676: {
677: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
678: PetscScalar *v;
680: PetscInt i;
683: *ncols = A->cmap->n;
684: if (cols) {
685: PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);
686: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
687: }
688: if (vals) {
689: PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);
690: v = mat->v + row;
691: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
692: }
693: return(0);
694: }
698: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
699: {
703: if (cols) {PetscFree(*cols);}
704: if (vals) {PetscFree(*vals); }
705: return(0);
706: }
707: /* ----------------------------------------------------------------*/
710: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
711: {
712: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
713: PetscInt i,j,idx=0;
717: if (!mat->roworiented) {
718: if (addv == INSERT_VALUES) {
719: for (j=0; j<n; j++) {
720: if (indexn[j] < 0) {idx += m; continue;}
721: #if defined(PETSC_USE_DEBUG)
722: 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);
723: #endif
724: for (i=0; i<m; i++) {
725: if (indexm[i] < 0) {idx++; continue;}
726: #if defined(PETSC_USE_DEBUG)
727: 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);
728: #endif
729: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
730: }
731: }
732: } else {
733: for (j=0; j<n; j++) {
734: if (indexn[j] < 0) {idx += m; continue;}
735: #if defined(PETSC_USE_DEBUG)
736: 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);
737: #endif
738: for (i=0; i<m; i++) {
739: if (indexm[i] < 0) {idx++; continue;}
740: #if defined(PETSC_USE_DEBUG)
741: 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);
742: #endif
743: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
744: }
745: }
746: }
747: } else {
748: if (addv == INSERT_VALUES) {
749: for (i=0; i<m; i++) {
750: if (indexm[i] < 0) { idx += n; continue;}
751: #if defined(PETSC_USE_DEBUG)
752: 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);
753: #endif
754: for (j=0; j<n; j++) {
755: if (indexn[j] < 0) { idx++; continue;}
756: #if defined(PETSC_USE_DEBUG)
757: 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);
758: #endif
759: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
760: }
761: }
762: } else {
763: for (i=0; i<m; i++) {
764: if (indexm[i] < 0) { idx += n; continue;}
765: #if defined(PETSC_USE_DEBUG)
766: 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);
767: #endif
768: for (j=0; j<n; j++) {
769: if (indexn[j] < 0) { idx++; continue;}
770: #if defined(PETSC_USE_DEBUG)
771: 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);
772: #endif
773: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
774: }
775: }
776: }
777: }
778: return(0);
779: }
783: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
784: {
785: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
786: PetscInt i,j;
789: /* row-oriented output */
790: for (i=0; i<m; i++) {
791: if (indexm[i] < 0) {v += n;continue;}
792: 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);
793: for (j=0; j<n; j++) {
794: if (indexn[j] < 0) {v++; continue;}
795: 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);
796: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
797: }
798: }
799: return(0);
800: }
802: /* -----------------------------------------------------------------*/
806: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
807: {
808: Mat_SeqDense *a;
810: PetscInt *scols,i,j,nz,header[4];
811: int fd;
812: PetscMPIInt size;
813: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
814: PetscScalar *vals,*svals,*v,*w;
815: MPI_Comm comm;
818: PetscObjectGetComm((PetscObject)viewer,&comm);
819: MPI_Comm_size(comm,&size);
820: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
821: PetscViewerBinaryGetDescriptor(viewer,&fd);
822: PetscBinaryRead(fd,header,4,PETSC_INT);
823: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
824: M = header[1]; N = header[2]; nz = header[3];
826: /* set global size if not set already*/
827: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
828: MatSetSizes(newmat,M,N,M,N);
829: } else {
830: /* if sizes and type are already set, check if the vector global sizes are correct */
831: MatGetSize(newmat,&grows,&gcols);
832: 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);
833: }
834: MatSeqDenseSetPreallocation(newmat,NULL);
836: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
837: a = (Mat_SeqDense*)newmat->data;
838: v = a->v;
839: /* Allocate some temp space to read in the values and then flip them
840: from row major to column major */
841: PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);
842: /* read in nonzero values */
843: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
844: /* now flip the values and store them in the matrix*/
845: for (j=0; j<N; j++) {
846: for (i=0; i<M; i++) {
847: *v++ =w[i*N+j];
848: }
849: }
850: PetscFree(w);
851: } else {
852: /* read row lengths */
853: PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);
854: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
856: a = (Mat_SeqDense*)newmat->data;
857: v = a->v;
859: /* read column indices and nonzeros */
860: PetscMalloc((nz+1)*sizeof(PetscInt),&scols);
861: cols = scols;
862: PetscBinaryRead(fd,cols,nz,PETSC_INT);
863: PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);
864: vals = svals;
865: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
867: /* insert into matrix */
868: for (i=0; i<M; i++) {
869: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
870: svals += rowlengths[i]; scols += rowlengths[i];
871: }
872: PetscFree(vals);
873: PetscFree(cols);
874: PetscFree(rowlengths);
875: }
876: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
877: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
878: return(0);
879: }
883: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
884: {
885: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
886: PetscErrorCode ierr;
887: PetscInt i,j;
888: const char *name;
889: PetscScalar *v;
890: PetscViewerFormat format;
891: #if defined(PETSC_USE_COMPLEX)
892: PetscBool allreal = PETSC_TRUE;
893: #endif
896: PetscViewerGetFormat(viewer,&format);
897: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
898: return(0); /* do nothing for now */
899: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
900: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
901: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
902: for (i=0; i<A->rmap->n; i++) {
903: v = a->v + i;
904: PetscViewerASCIIPrintf(viewer,"row %D:",i);
905: for (j=0; j<A->cmap->n; j++) {
906: #if defined(PETSC_USE_COMPLEX)
907: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
908: PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));
909: } else if (PetscRealPart(*v)) {
910: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));
911: }
912: #else
913: if (*v) {
914: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);
915: }
916: #endif
917: v += a->lda;
918: }
919: PetscViewerASCIIPrintf(viewer,"\n");
920: }
921: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
922: } else {
923: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
924: #if defined(PETSC_USE_COMPLEX)
925: /* determine if matrix has all real values */
926: v = a->v;
927: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
928: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
929: }
930: #endif
931: if (format == PETSC_VIEWER_ASCII_MATLAB) {
932: PetscObjectGetName((PetscObject)A,&name);
933: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
934: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
935: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
936: } else {
937: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
938: }
940: for (i=0; i<A->rmap->n; i++) {
941: v = a->v + i;
942: for (j=0; j<A->cmap->n; j++) {
943: #if defined(PETSC_USE_COMPLEX)
944: if (allreal) {
945: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
946: } else {
947: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
948: }
949: #else
950: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
951: #endif
952: v += a->lda;
953: }
954: PetscViewerASCIIPrintf(viewer,"\n");
955: }
956: if (format == PETSC_VIEWER_ASCII_MATLAB) {
957: PetscViewerASCIIPrintf(viewer,"];\n");
958: }
959: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
960: }
961: PetscViewerFlush(viewer);
962: return(0);
963: }
967: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
968: {
969: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
970: PetscErrorCode ierr;
971: int fd;
972: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
973: PetscScalar *v,*anonz,*vals;
974: PetscViewerFormat format;
977: PetscViewerBinaryGetDescriptor(viewer,&fd);
979: PetscViewerGetFormat(viewer,&format);
980: if (format == PETSC_VIEWER_NATIVE) {
981: /* store the matrix as a dense matrix */
982: PetscMalloc(4*sizeof(PetscInt),&col_lens);
984: col_lens[0] = MAT_FILE_CLASSID;
985: col_lens[1] = m;
986: col_lens[2] = n;
987: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
989: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
990: PetscFree(col_lens);
992: /* write out matrix, by rows */
993: PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);
994: v = a->v;
995: for (j=0; j<n; j++) {
996: for (i=0; i<m; i++) {
997: vals[j + i*n] = *v++;
998: }
999: }
1000: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1001: PetscFree(vals);
1002: } else {
1003: PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);
1005: col_lens[0] = MAT_FILE_CLASSID;
1006: col_lens[1] = m;
1007: col_lens[2] = n;
1008: col_lens[3] = nz;
1010: /* store lengths of each row and write (including header) to file */
1011: for (i=0; i<m; i++) col_lens[4+i] = n;
1012: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1014: /* Possibly should write in smaller increments, not whole matrix at once? */
1015: /* store column indices (zero start index) */
1016: ict = 0;
1017: for (i=0; i<m; i++) {
1018: for (j=0; j<n; j++) col_lens[ict++] = j;
1019: }
1020: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1021: PetscFree(col_lens);
1023: /* store nonzero values */
1024: PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);
1025: ict = 0;
1026: for (i=0; i<m; i++) {
1027: v = a->v + i;
1028: for (j=0; j<n; j++) {
1029: anonz[ict++] = *v; v += a->lda;
1030: }
1031: }
1032: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1033: PetscFree(anonz);
1034: }
1035: return(0);
1036: }
1038: #include <petscdraw.h>
1041: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1042: {
1043: Mat A = (Mat) Aa;
1044: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1045: PetscErrorCode ierr;
1046: PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j;
1047: PetscScalar *v = a->v;
1048: PetscViewer viewer;
1049: PetscDraw popup;
1050: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1051: PetscViewerFormat format;
1054: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1055: PetscViewerGetFormat(viewer,&format);
1056: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1058: /* Loop over matrix elements drawing boxes */
1059: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1060: /* Blue for negative and Red for positive */
1061: color = PETSC_DRAW_BLUE;
1062: for (j = 0; j < n; j++) {
1063: x_l = j;
1064: x_r = x_l + 1.0;
1065: for (i = 0; i < m; i++) {
1066: y_l = m - i - 1.0;
1067: y_r = y_l + 1.0;
1068: if (PetscRealPart(v[j*m+i]) > 0.) {
1069: color = PETSC_DRAW_RED;
1070: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1071: color = PETSC_DRAW_BLUE;
1072: } else {
1073: continue;
1074: }
1075: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1076: }
1077: }
1078: } else {
1079: /* use contour shading to indicate magnitude of values */
1080: /* first determine max of all nonzero values */
1081: for (i = 0; i < m*n; i++) {
1082: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1083: }
1084: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1085: PetscDrawGetPopup(draw,&popup);
1086: if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1087: for (j = 0; j < n; j++) {
1088: x_l = j;
1089: x_r = x_l + 1.0;
1090: for (i = 0; i < m; i++) {
1091: y_l = m - i - 1.0;
1092: y_r = y_l + 1.0;
1093: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1094: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1095: }
1096: }
1097: }
1098: return(0);
1099: }
1103: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1104: {
1105: PetscDraw draw;
1106: PetscBool isnull;
1107: PetscReal xr,yr,xl,yl,h,w;
1111: PetscViewerDrawGetDraw(viewer,0,&draw);
1112: PetscDrawIsNull(draw,&isnull);
1113: if (isnull) return(0);
1115: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1116: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1117: xr += w; yr += h; xl = -w; yl = -h;
1118: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1119: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1120: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1121: return(0);
1122: }
1126: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1127: {
1129: PetscBool iascii,isbinary,isdraw;
1132: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1133: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1134: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1136: if (iascii) {
1137: MatView_SeqDense_ASCII(A,viewer);
1138: } else if (isbinary) {
1139: MatView_SeqDense_Binary(A,viewer);
1140: } else if (isdraw) {
1141: MatView_SeqDense_Draw(A,viewer);
1142: }
1143: return(0);
1144: }
1148: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1149: {
1150: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1154: #if defined(PETSC_USE_LOG)
1155: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1156: #endif
1157: PetscFree(l->pivots);
1158: if (!l->user_alloc) {PetscFree(l->v);}
1159: PetscFree(mat->data);
1161: PetscObjectChangeTypeName((PetscObject)mat,0);
1162: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1163: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1164: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1165: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1166: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1167: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1168: return(0);
1169: }
1173: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1174: {
1175: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1177: PetscInt k,j,m,n,M;
1178: PetscScalar *v,tmp;
1181: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1182: if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1183: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1184: else {
1185: for (j=0; j<m; j++) {
1186: for (k=0; k<j; k++) {
1187: tmp = v[j + k*M];
1188: v[j + k*M] = v[k + j*M];
1189: v[k + j*M] = tmp;
1190: }
1191: }
1192: }
1193: } else { /* out-of-place transpose */
1194: Mat tmat;
1195: Mat_SeqDense *tmatd;
1196: PetscScalar *v2;
1198: if (reuse == MAT_INITIAL_MATRIX) {
1199: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1200: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1201: MatSetType(tmat,((PetscObject)A)->type_name);
1202: MatSeqDenseSetPreallocation(tmat,NULL);
1203: } else {
1204: tmat = *matout;
1205: }
1206: tmatd = (Mat_SeqDense*)tmat->data;
1207: v = mat->v; v2 = tmatd->v;
1208: for (j=0; j<n; j++) {
1209: for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
1210: }
1211: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1212: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1213: *matout = tmat;
1214: }
1215: return(0);
1216: }
1220: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1221: {
1222: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1223: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1224: PetscInt i,j;
1225: PetscScalar *v1,*v2;
1228: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1229: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1230: for (i=0; i<A1->rmap->n; i++) {
1231: v1 = mat1->v+i; v2 = mat2->v+i;
1232: for (j=0; j<A1->cmap->n; j++) {
1233: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1234: v1 += mat1->lda; v2 += mat2->lda;
1235: }
1236: }
1237: *flg = PETSC_TRUE;
1238: return(0);
1239: }
1243: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1244: {
1245: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1247: PetscInt i,n,len;
1248: PetscScalar *x,zero = 0.0;
1251: VecSet(v,zero);
1252: VecGetSize(v,&n);
1253: VecGetArray(v,&x);
1254: len = PetscMin(A->rmap->n,A->cmap->n);
1255: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1256: for (i=0; i<len; i++) {
1257: x[i] = mat->v[i*mat->lda + i];
1258: }
1259: VecRestoreArray(v,&x);
1260: return(0);
1261: }
1265: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1266: {
1267: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1268: PetscScalar *l,*r,x,*v;
1270: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1273: if (ll) {
1274: VecGetSize(ll,&m);
1275: VecGetArray(ll,&l);
1276: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1277: for (i=0; i<m; i++) {
1278: x = l[i];
1279: v = mat->v + i;
1280: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1281: }
1282: VecRestoreArray(ll,&l);
1283: PetscLogFlops(n*m);
1284: }
1285: if (rr) {
1286: VecGetSize(rr,&n);
1287: VecGetArray(rr,&r);
1288: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1289: for (i=0; i<n; i++) {
1290: x = r[i];
1291: v = mat->v + i*m;
1292: for (j=0; j<m; j++) (*v++) *= x;
1293: }
1294: VecRestoreArray(rr,&r);
1295: PetscLogFlops(n*m);
1296: }
1297: return(0);
1298: }
1302: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1303: {
1304: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1305: PetscScalar *v = mat->v;
1306: PetscReal sum = 0.0;
1307: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1311: if (type == NORM_FROBENIUS) {
1312: if (lda>m) {
1313: for (j=0; j<A->cmap->n; j++) {
1314: v = mat->v+j*lda;
1315: for (i=0; i<m; i++) {
1316: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1317: }
1318: }
1319: } else {
1320: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1321: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1322: }
1323: }
1324: *nrm = PetscSqrtReal(sum);
1325: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1326: } else if (type == NORM_1) {
1327: *nrm = 0.0;
1328: for (j=0; j<A->cmap->n; j++) {
1329: v = mat->v + j*mat->lda;
1330: sum = 0.0;
1331: for (i=0; i<A->rmap->n; i++) {
1332: sum += PetscAbsScalar(*v); v++;
1333: }
1334: if (sum > *nrm) *nrm = sum;
1335: }
1336: PetscLogFlops(A->cmap->n*A->rmap->n);
1337: } else if (type == NORM_INFINITY) {
1338: *nrm = 0.0;
1339: for (j=0; j<A->rmap->n; j++) {
1340: v = mat->v + j;
1341: sum = 0.0;
1342: for (i=0; i<A->cmap->n; i++) {
1343: sum += PetscAbsScalar(*v); v += mat->lda;
1344: }
1345: if (sum > *nrm) *nrm = sum;
1346: }
1347: PetscLogFlops(A->cmap->n*A->rmap->n);
1348: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1349: return(0);
1350: }
1354: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1355: {
1356: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1360: switch (op) {
1361: case MAT_ROW_ORIENTED:
1362: aij->roworiented = flg;
1363: break;
1364: case MAT_NEW_NONZERO_LOCATIONS:
1365: case MAT_NEW_NONZERO_LOCATION_ERR:
1366: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1367: case MAT_NEW_DIAGONALS:
1368: case MAT_KEEP_NONZERO_PATTERN:
1369: case MAT_IGNORE_OFF_PROC_ENTRIES:
1370: case MAT_USE_HASH_TABLE:
1371: case MAT_IGNORE_LOWER_TRIANGULAR:
1372: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1373: break;
1374: case MAT_SPD:
1375: case MAT_SYMMETRIC:
1376: case MAT_STRUCTURALLY_SYMMETRIC:
1377: case MAT_HERMITIAN:
1378: case MAT_SYMMETRY_ETERNAL:
1379: /* These options are handled directly by MatSetOption() */
1380: break;
1381: default:
1382: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1383: }
1384: return(0);
1385: }
1389: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1390: {
1391: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1393: PetscInt lda=l->lda,m=A->rmap->n,j;
1396: if (lda>m) {
1397: for (j=0; j<A->cmap->n; j++) {
1398: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1399: }
1400: } else {
1401: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1402: }
1403: return(0);
1404: }
1408: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1409: {
1410: PetscErrorCode ierr;
1411: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1412: PetscInt m = l->lda, n = A->cmap->n, i,j;
1413: PetscScalar *slot,*bb;
1414: const PetscScalar *xx;
1417: #if defined(PETSC_USE_DEBUG)
1418: for (i=0; i<N; i++) {
1419: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1420: 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);
1421: }
1422: #endif
1424: /* fix right hand side if needed */
1425: if (x && b) {
1426: VecGetArrayRead(x,&xx);
1427: VecGetArray(b,&bb);
1428: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1429: VecRestoreArrayRead(x,&xx);
1430: VecRestoreArray(b,&bb);
1431: }
1433: for (i=0; i<N; i++) {
1434: slot = l->v + rows[i];
1435: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1436: }
1437: if (diag != 0.0) {
1438: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1439: for (i=0; i<N; i++) {
1440: slot = l->v + (m+1)*rows[i];
1441: *slot = diag;
1442: }
1443: }
1444: return(0);
1445: }
1449: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1450: {
1451: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1454: 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");
1455: *array = mat->v;
1456: return(0);
1457: }
1461: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1462: {
1464: *array = 0; /* user cannot accidently use the array later */
1465: return(0);
1466: }
1470: /*@C
1471: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1473: Not Collective
1475: Input Parameter:
1476: . mat - a MATSEQDENSE matrix
1478: Output Parameter:
1479: . array - pointer to the data
1481: Level: intermediate
1483: .seealso: MatDenseRestoreArray()
1484: @*/
1485: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1486: {
1490: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1491: return(0);
1492: }
1496: /*@C
1497: MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray()
1499: Not Collective
1501: Input Parameters:
1502: . mat - a MATSEQDENSE matrix
1503: . array - pointer to the data
1505: Level: intermediate
1507: .seealso: MatDenseGetArray()
1508: @*/
1509: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1510: {
1514: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1515: return(0);
1516: }
1520: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1521: {
1522: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1524: PetscInt i,j,nrows,ncols;
1525: const PetscInt *irow,*icol;
1526: PetscScalar *av,*bv,*v = mat->v;
1527: Mat newmat;
1530: ISGetIndices(isrow,&irow);
1531: ISGetIndices(iscol,&icol);
1532: ISGetLocalSize(isrow,&nrows);
1533: ISGetLocalSize(iscol,&ncols);
1535: /* Check submatrixcall */
1536: if (scall == MAT_REUSE_MATRIX) {
1537: PetscInt n_cols,n_rows;
1538: MatGetSize(*B,&n_rows,&n_cols);
1539: if (n_rows != nrows || n_cols != ncols) {
1540: /* resize the result matrix to match number of requested rows/columns */
1541: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1542: }
1543: newmat = *B;
1544: } else {
1545: /* Create and fill new matrix */
1546: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1547: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1548: MatSetType(newmat,((PetscObject)A)->type_name);
1549: MatSeqDenseSetPreallocation(newmat,NULL);
1550: }
1552: /* Now extract the data pointers and do the copy,column at a time */
1553: bv = ((Mat_SeqDense*)newmat->data)->v;
1555: for (i=0; i<ncols; i++) {
1556: av = v + mat->lda*icol[i];
1557: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1558: }
1560: /* Assemble the matrices so that the correct flags are set */
1561: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1562: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1564: /* Free work space */
1565: ISRestoreIndices(isrow,&irow);
1566: ISRestoreIndices(iscol,&icol);
1567: *B = newmat;
1568: return(0);
1569: }
1573: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1574: {
1576: PetscInt i;
1579: if (scall == MAT_INITIAL_MATRIX) {
1580: PetscMalloc((n+1)*sizeof(Mat),B);
1581: }
1583: for (i=0; i<n; i++) {
1584: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1585: }
1586: return(0);
1587: }
1591: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1592: {
1594: return(0);
1595: }
1599: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1600: {
1602: return(0);
1603: }
1607: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1608: {
1609: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1611: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1614: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1615: if (A->ops->copy != B->ops->copy) {
1616: MatCopy_Basic(A,B,str);
1617: return(0);
1618: }
1619: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1620: if (lda1>m || lda2>m) {
1621: for (j=0; j<n; j++) {
1622: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1623: }
1624: } else {
1625: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1626: }
1627: return(0);
1628: }
1632: PetscErrorCode MatSetUp_SeqDense(Mat A)
1633: {
1637: MatSeqDenseSetPreallocation(A,0);
1638: return(0);
1639: }
1643: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1644: {
1645: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1646: PetscInt i,nz = A->rmap->n*A->cmap->n;
1647: PetscScalar *aa = a->v;
1650: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1651: return(0);
1652: }
1656: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1657: {
1658: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1659: PetscInt i,nz = A->rmap->n*A->cmap->n;
1660: PetscScalar *aa = a->v;
1663: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1664: return(0);
1665: }
1669: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1670: {
1671: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1672: PetscInt i,nz = A->rmap->n*A->cmap->n;
1673: PetscScalar *aa = a->v;
1676: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1677: return(0);
1678: }
1680: /* ----------------------------------------------------------------*/
1683: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1684: {
1688: if (scall == MAT_INITIAL_MATRIX) {
1689: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1690: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1691: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1692: }
1693: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1694: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1695: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1696: return(0);
1697: }
1701: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1702: {
1704: PetscInt m=A->rmap->n,n=B->cmap->n;
1705: Mat Cmat;
1708: 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);
1709: MatCreate(PETSC_COMM_SELF,&Cmat);
1710: MatSetSizes(Cmat,m,n,m,n);
1711: MatSetType(Cmat,MATSEQDENSE);
1712: MatSeqDenseSetPreallocation(Cmat,NULL);
1714: *C = Cmat;
1715: return(0);
1716: }
1720: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1721: {
1722: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1723: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1724: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1725: PetscBLASInt m,n,k;
1726: PetscScalar _DOne=1.0,_DZero=0.0;
1730: PetscBLASIntCast(A->rmap->n,&m);
1731: PetscBLASIntCast(B->cmap->n,&n);
1732: PetscBLASIntCast(A->cmap->n,&k);
1733: PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1734: return(0);
1735: }
1739: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1740: {
1744: if (scall == MAT_INITIAL_MATRIX) {
1745: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
1746: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1747: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
1748: }
1749: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
1750: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1751: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
1752: return(0);
1753: }
1757: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1758: {
1760: PetscInt m=A->cmap->n,n=B->cmap->n;
1761: Mat Cmat;
1764: 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);
1765: MatCreate(PETSC_COMM_SELF,&Cmat);
1766: MatSetSizes(Cmat,m,n,m,n);
1767: MatSetType(Cmat,MATSEQDENSE);
1768: MatSeqDenseSetPreallocation(Cmat,NULL);
1770: Cmat->assembled = PETSC_TRUE;
1772: *C = Cmat;
1773: return(0);
1774: }
1778: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1779: {
1780: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1781: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1782: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1783: PetscBLASInt m,n,k;
1784: PetscScalar _DOne=1.0,_DZero=0.0;
1788: PetscBLASIntCast(A->cmap->n,&m);
1789: PetscBLASIntCast(B->cmap->n,&n);
1790: PetscBLASIntCast(A->rmap->n,&k);
1791: /*
1792: Note the m and n arguments below are the number rows and columns of A', not A!
1793: */
1794: PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1795: return(0);
1796: }
1800: PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1801: {
1802: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1804: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1805: PetscScalar *x;
1806: MatScalar *aa = a->v;
1809: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1811: VecSet(v,0.0);
1812: VecGetArray(v,&x);
1813: VecGetLocalSize(v,&p);
1814: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1815: for (i=0; i<m; i++) {
1816: x[i] = aa[i]; if (idx) idx[i] = 0;
1817: for (j=1; j<n; j++) {
1818: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1819: }
1820: }
1821: VecRestoreArray(v,&x);
1822: return(0);
1823: }
1827: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1828: {
1829: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1831: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1832: PetscScalar *x;
1833: PetscReal atmp;
1834: MatScalar *aa = a->v;
1837: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1839: VecSet(v,0.0);
1840: VecGetArray(v,&x);
1841: VecGetLocalSize(v,&p);
1842: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1843: for (i=0; i<m; i++) {
1844: x[i] = PetscAbsScalar(aa[i]);
1845: for (j=1; j<n; j++) {
1846: atmp = PetscAbsScalar(aa[i+m*j]);
1847: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1848: }
1849: }
1850: VecRestoreArray(v,&x);
1851: return(0);
1852: }
1856: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1857: {
1858: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1860: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1861: PetscScalar *x;
1862: MatScalar *aa = a->v;
1865: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1867: VecSet(v,0.0);
1868: VecGetArray(v,&x);
1869: VecGetLocalSize(v,&p);
1870: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1871: for (i=0; i<m; i++) {
1872: x[i] = aa[i]; if (idx) idx[i] = 0;
1873: for (j=1; j<n; j++) {
1874: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1875: }
1876: }
1877: VecRestoreArray(v,&x);
1878: return(0);
1879: }
1883: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1884: {
1885: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1887: PetscScalar *x;
1890: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1892: VecGetArray(v,&x);
1893: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1894: VecRestoreArray(v,&x);
1895: return(0);
1896: }
1901: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1902: {
1904: PetscInt i,j,m,n;
1905: PetscScalar *a;
1908: MatGetSize(A,&m,&n);
1909: PetscMemzero(norms,n*sizeof(PetscReal));
1910: MatDenseGetArray(A,&a);
1911: if (type == NORM_2) {
1912: for (i=0; i<n; i++) {
1913: for (j=0; j<m; j++) {
1914: norms[i] += PetscAbsScalar(a[j]*a[j]);
1915: }
1916: a += m;
1917: }
1918: } else if (type == NORM_1) {
1919: for (i=0; i<n; i++) {
1920: for (j=0; j<m; j++) {
1921: norms[i] += PetscAbsScalar(a[j]);
1922: }
1923: a += m;
1924: }
1925: } else if (type == NORM_INFINITY) {
1926: for (i=0; i<n; i++) {
1927: for (j=0; j<m; j++) {
1928: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1929: }
1930: a += m;
1931: }
1932: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
1933: MatDenseRestoreArray(A,&a);
1934: if (type == NORM_2) {
1935: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1936: }
1937: return(0);
1938: }
1942: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
1943: {
1945: PetscScalar *a;
1946: PetscInt m,n,i;
1949: MatGetSize(x,&m,&n);
1950: MatDenseGetArray(x,&a);
1951: for (i=0; i<m*n; i++) {
1952: PetscRandomGetValue(rctx,a+i);
1953: }
1954: MatDenseRestoreArray(x,&a);
1955: return(0);
1956: }
1959: /* -------------------------------------------------------------------*/
1960: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
1961: MatGetRow_SeqDense,
1962: MatRestoreRow_SeqDense,
1963: MatMult_SeqDense,
1964: /* 4*/ MatMultAdd_SeqDense,
1965: MatMultTranspose_SeqDense,
1966: MatMultTransposeAdd_SeqDense,
1967: 0,
1968: 0,
1969: 0,
1970: /* 10*/ 0,
1971: MatLUFactor_SeqDense,
1972: MatCholeskyFactor_SeqDense,
1973: MatSOR_SeqDense,
1974: MatTranspose_SeqDense,
1975: /* 15*/ MatGetInfo_SeqDense,
1976: MatEqual_SeqDense,
1977: MatGetDiagonal_SeqDense,
1978: MatDiagonalScale_SeqDense,
1979: MatNorm_SeqDense,
1980: /* 20*/ MatAssemblyBegin_SeqDense,
1981: MatAssemblyEnd_SeqDense,
1982: MatSetOption_SeqDense,
1983: MatZeroEntries_SeqDense,
1984: /* 24*/ MatZeroRows_SeqDense,
1985: 0,
1986: 0,
1987: 0,
1988: 0,
1989: /* 29*/ MatSetUp_SeqDense,
1990: 0,
1991: 0,
1992: 0,
1993: 0,
1994: /* 34*/ MatDuplicate_SeqDense,
1995: 0,
1996: 0,
1997: 0,
1998: 0,
1999: /* 39*/ MatAXPY_SeqDense,
2000: MatGetSubMatrices_SeqDense,
2001: 0,
2002: MatGetValues_SeqDense,
2003: MatCopy_SeqDense,
2004: /* 44*/ MatGetRowMax_SeqDense,
2005: MatScale_SeqDense,
2006: 0,
2007: 0,
2008: 0,
2009: /* 49*/ MatSetRandom_SeqDense,
2010: 0,
2011: 0,
2012: 0,
2013: 0,
2014: /* 54*/ 0,
2015: 0,
2016: 0,
2017: 0,
2018: 0,
2019: /* 59*/ 0,
2020: MatDestroy_SeqDense,
2021: MatView_SeqDense,
2022: 0,
2023: 0,
2024: /* 64*/ 0,
2025: 0,
2026: 0,
2027: 0,
2028: 0,
2029: /* 69*/ MatGetRowMaxAbs_SeqDense,
2030: 0,
2031: 0,
2032: 0,
2033: 0,
2034: /* 74*/ 0,
2035: 0,
2036: 0,
2037: 0,
2038: 0,
2039: /* 79*/ 0,
2040: 0,
2041: 0,
2042: 0,
2043: /* 83*/ MatLoad_SeqDense,
2044: 0,
2045: MatIsHermitian_SeqDense,
2046: 0,
2047: 0,
2048: 0,
2049: /* 89*/ MatMatMult_SeqDense_SeqDense,
2050: MatMatMultSymbolic_SeqDense_SeqDense,
2051: MatMatMultNumeric_SeqDense_SeqDense,
2052: 0,
2053: 0,
2054: /* 94*/ 0,
2055: 0,
2056: 0,
2057: 0,
2058: 0,
2059: /* 99*/ 0,
2060: 0,
2061: 0,
2062: MatConjugate_SeqDense,
2063: 0,
2064: /*104*/ 0,
2065: MatRealPart_SeqDense,
2066: MatImaginaryPart_SeqDense,
2067: 0,
2068: 0,
2069: /*109*/ MatMatSolve_SeqDense,
2070: 0,
2071: MatGetRowMin_SeqDense,
2072: MatGetColumnVector_SeqDense,
2073: 0,
2074: /*114*/ 0,
2075: 0,
2076: 0,
2077: 0,
2078: 0,
2079: /*119*/ 0,
2080: 0,
2081: 0,
2082: 0,
2083: 0,
2084: /*124*/ 0,
2085: MatGetColumnNorms_SeqDense,
2086: 0,
2087: 0,
2088: 0,
2089: /*129*/ 0,
2090: MatTransposeMatMult_SeqDense_SeqDense,
2091: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2092: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2093: 0,
2094: /*134*/ 0,
2095: 0,
2096: 0,
2097: 0,
2098: 0,
2099: /*139*/ 0,
2100: 0
2101: };
2105: /*@C
2106: MatCreateSeqDense - Creates a sequential dense matrix that
2107: is stored in column major order (the usual Fortran 77 manner). Many
2108: of the matrix operations use the BLAS and LAPACK routines.
2110: Collective on MPI_Comm
2112: Input Parameters:
2113: + comm - MPI communicator, set to PETSC_COMM_SELF
2114: . m - number of rows
2115: . n - number of columns
2116: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2117: to control all matrix memory allocation.
2119: Output Parameter:
2120: . A - the matrix
2122: Notes:
2123: The data input variable is intended primarily for Fortran programmers
2124: who wish to allocate their own matrix memory space. Most users should
2125: set data=NULL.
2127: Level: intermediate
2129: .keywords: dense, matrix, LAPACK, BLAS
2131: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2132: @*/
2133: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2134: {
2138: MatCreate(comm,A);
2139: MatSetSizes(*A,m,n,m,n);
2140: MatSetType(*A,MATSEQDENSE);
2141: MatSeqDenseSetPreallocation(*A,data);
2142: return(0);
2143: }
2147: /*@C
2148: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2150: Collective on MPI_Comm
2152: Input Parameters:
2153: + A - the matrix
2154: - data - the array (or NULL)
2156: Notes:
2157: The data input variable is intended primarily for Fortran programmers
2158: who wish to allocate their own matrix memory space. Most users should
2159: need not call this routine.
2161: Level: intermediate
2163: .keywords: dense, matrix, LAPACK, BLAS
2165: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2167: @*/
2168: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2169: {
2173: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2174: return(0);
2175: }
2179: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2180: {
2181: Mat_SeqDense *b;
2185: B->preallocated = PETSC_TRUE;
2187: PetscLayoutSetUp(B->rmap);
2188: PetscLayoutSetUp(B->cmap);
2190: b = (Mat_SeqDense*)B->data;
2191: b->Mmax = B->rmap->n;
2192: b->Nmax = B->cmap->n;
2193: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2195: if (!data) { /* petsc-allocated storage */
2196: if (!b->user_alloc) { PetscFree(b->v); }
2197: PetscMalloc((size_t)b->lda*b->Nmax*sizeof(PetscScalar),&b->v);
2198: PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));
2199: PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));
2201: b->user_alloc = PETSC_FALSE;
2202: } else { /* user-allocated storage */
2203: if (!b->user_alloc) { PetscFree(b->v); }
2204: b->v = data;
2205: b->user_alloc = PETSC_TRUE;
2206: }
2207: B->assembled = PETSC_TRUE;
2208: return(0);
2209: }
2213: /*@C
2214: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2216: Input parameter:
2217: + A - the matrix
2218: - lda - the leading dimension
2220: Notes:
2221: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2222: it asserts that the preallocation has a leading dimension (the LDA parameter
2223: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2225: Level: intermediate
2227: .keywords: dense, matrix, LAPACK, BLAS
2229: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2231: @*/
2232: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2233: {
2234: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2237: 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);
2238: b->lda = lda;
2239: b->changelda = PETSC_FALSE;
2240: b->Mmax = PetscMax(b->Mmax,lda);
2241: return(0);
2242: }
2244: /*MC
2245: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2247: Options Database Keys:
2248: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2250: Level: beginner
2252: .seealso: MatCreateSeqDense()
2254: M*/
2258: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2259: {
2260: Mat_SeqDense *b;
2262: PetscMPIInt size;
2265: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2266: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2268: PetscNewLog(B,Mat_SeqDense,&b);
2269: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2270: B->data = (void*)b;
2272: b->pivots = 0;
2273: b->roworiented = PETSC_TRUE;
2274: b->v = 0;
2275: b->changelda = PETSC_FALSE;
2277: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2278: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2280: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2281: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);
2282: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2283: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2284: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2285: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2286: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2287: return(0);
2288: }