Actual source code: dense.c
petsc-3.3-p7 2013-05-11
2: /*
3: Defines the basic matrix operations for sequential dense.
4: */
6: #include <../src/mat/impls/dense/seq/dense.h>
7: #include <petscblaslapack.h>
9: #include <../src/mat/impls/aij/seq/aij.h>
10: EXTERN_C_BEGIN
13: 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(((PetscObject)A)->comm,&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,PETSC_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);
38:
39: if (reuse == MAT_REUSE_MATRIX) {
40: MatHeaderReplace(A,B);
41: } else {
42: *newmat = B;
43: }
44: return(0);
45: }
46: EXTERN_C_END
50: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
51: {
52: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
53: PetscScalar oalpha = alpha;
54: PetscInt j;
55: PetscBLASInt N,m,ldax,lday,one = 1;
59: N = PetscBLASIntCast(X->rmap->n*X->cmap->n);
60: m = PetscBLASIntCast(X->rmap->n);
61: ldax = PetscBLASIntCast(x->lda);
62: lday = PetscBLASIntCast(y->lda);
63: if (ldax>m || lday>m) {
64: for (j=0; j<X->cmap->n; j++) {
65: BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one);
66: }
67: } else {
68: BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one);
69: }
70: PetscLogFlops(PetscMax(2*N-1,0));
71: return(0);
72: }
76: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
77: {
78: PetscInt N = A->rmap->n*A->cmap->n;
81: info->block_size = 1.0;
82: info->nz_allocated = (double)N;
83: info->nz_used = (double)N;
84: info->nz_unneeded = (double)0;
85: info->assemblies = (double)A->num_ass;
86: info->mallocs = 0;
87: info->memory = ((PetscObject)A)->mem;
88: info->fill_ratio_given = 0;
89: info->fill_ratio_needed = 0;
90: info->factor_mallocs = 0;
91: return(0);
92: }
96: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
97: {
98: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
99: PetscScalar oalpha = alpha;
101: PetscBLASInt one = 1,j,nz,lda = PetscBLASIntCast(a->lda);
104: if (lda>A->rmap->n) {
105: nz = PetscBLASIntCast(A->rmap->n);
106: for (j=0; j<A->cmap->n; j++) {
107: BLASscal_(&nz,&oalpha,a->v+j*lda,&one);
108: }
109: } else {
110: nz = PetscBLASIntCast(A->rmap->n*A->cmap->n);
111: 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: }
138:
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,PETSC_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(((PetscObject)A)->comm,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 = PetscBLASIntCast(A->rmap->n);
204:
206: VecGetArray(xx,&x);
207: VecGetArray(yy,&y);
208: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
209: if (A->factortype == MAT_FACTOR_LU) {
210: #if defined(PETSC_MISSING_LAPACK_GETRS)
211: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
212: #else
213: LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
214: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
215: #endif
216: } else if (A->factortype == MAT_FACTOR_CHOLESKY){
217: #if defined(PETSC_MISSING_LAPACK_POTRS)
218: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
219: #else
220: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
221: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
222: #endif
223: }
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: m = PetscBLASIntCast(A->rmap->n);
244: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
245: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
246: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
247: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
249: MatGetSize(B,PETSC_NULL,&n);
250: nrhs = PetscBLASIntCast(n);
251: MatGetArray(B,&b);
252: MatGetArray(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: 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: 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: }
271: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
273: MatRestoreArray(B,&b);
274: MatRestoreArray(X,&x);
275: PetscLogFlops(nrhs*(2.0*m*m - m));
276: return(0);
277: }
281: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
282: {
283: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
285: PetscScalar *x,*y;
286: PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n);
287:
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: 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: 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 = PetscBLASIntCast(A->rmap->n);
323:
325: VecGetArray(xx,&x);
326: VecGetArray(yy,&y);
327: if (!A->rmap->n || !A->cmap->n) return(0);
328: if (yy == zz) {
329: VecDuplicate(yy,&tmp);
330: PetscLogObjectParent(A,tmp);
331: VecCopy(yy,tmp);
332: }
333: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
334: /* assume if pivots exist then use LU; else Cholesky */
335: if (mat->pivots) {
336: #if defined(PETSC_MISSING_LAPACK_GETRS)
337: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
338: #else
339: LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
340: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
341: #endif
342: } else {
343: #if defined(PETSC_MISSING_LAPACK_POTRS)
344: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
345: #else
346: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
347: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
348: #endif
349: }
350: if (tmp) {
351: VecAXPY(yy,sone,tmp);
352: VecDestroy(&tmp);
353: } else {
354: VecAXPY(yy,sone,zz);
355: }
356: VecRestoreArray(xx,&x);
357: VecRestoreArray(yy,&y);
358: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
359: return(0);
360: }
364: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
365: {
366: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
368: PetscScalar *x,*y,sone = 1.0;
369: Vec tmp;
370: PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n);
371:
373: if (!A->rmap->n || !A->cmap->n) return(0);
374: VecGetArray(xx,&x);
375: VecGetArray(yy,&y);
376: if (yy == zz) {
377: VecDuplicate(yy,&tmp);
378: PetscLogObjectParent(A,tmp);
379: VecCopy(yy,tmp);
380: }
381: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
382: /* assume if pivots exist then use LU; else Cholesky */
383: if (mat->pivots) {
384: #if defined(PETSC_MISSING_LAPACK_GETRS)
385: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
386: #else
387: LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
388: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
389: #endif
390: } else {
391: #if defined(PETSC_MISSING_LAPACK_POTRS)
392: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
393: #else
394: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
395: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve");
396: #endif
397: }
398: if (tmp) {
399: VecAXPY(yy,sone,tmp);
400: VecDestroy(&tmp);
401: } else {
402: VecAXPY(yy,sone,zz);
403: }
404: VecRestoreArray(xx,&x);
405: VecRestoreArray(yy,&y);
406: PetscLogFlops(2.0*A->cmap->n*A->cmap->n);
407: return(0);
408: }
410: /* ---------------------------------------------------------------*/
411: /* COMMENT: I have chosen to hide row permutation in the pivots,
412: rather than put it in the Mat->row slot.*/
415: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
416: {
417: #if defined(PETSC_MISSING_LAPACK_GETRF)
419: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
420: #else
421: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
423: PetscBLASInt n,m,info;
426: n = PetscBLASIntCast(A->cmap->n);
427: m = PetscBLASIntCast(A->rmap->n);
428: if (!mat->pivots) {
429: PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);
430: PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));
431: }
432: if (!A->rmap->n || !A->cmap->n) return(0);
433: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
434: LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info);
435: PetscFPTrapPop();
437: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
438: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
439: A->ops->solve = MatSolve_SeqDense;
440: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
441: A->ops->solveadd = MatSolveAdd_SeqDense;
442: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
443: A->factortype = MAT_FACTOR_LU;
445: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
446: #endif
447: return(0);
448: }
452: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
453: {
454: #if defined(PETSC_MISSING_LAPACK_POTRF)
456: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
457: #else
458: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
460: PetscBLASInt info,n = PetscBLASIntCast(A->cmap->n);
461:
463: PetscFree(mat->pivots);
465: if (!A->rmap->n || !A->cmap->n) return(0);
466: LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info);
467: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
468: A->ops->solve = MatSolve_SeqDense;
469: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
470: A->ops->solveadd = MatSolveAdd_SeqDense;
471: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
472: A->factortype = MAT_FACTOR_CHOLESKY;
473: PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
474: #endif
475: return(0);
476: }
481: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
482: {
484: MatFactorInfo info;
487: info.fill = 1.0;
488: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
489: MatCholeskyFactor_SeqDense(fact,0,&info);
490: return(0);
491: }
495: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
496: {
498: fact->assembled = PETSC_TRUE;
499: fact->preallocated = PETSC_TRUE;
500: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
501: return(0);
502: }
506: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
507: {
509: fact->preallocated = PETSC_TRUE;
510: fact->assembled = PETSC_TRUE;
511: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
512: return(0);
513: }
515: EXTERN_C_BEGIN
518: PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
519: {
523: MatCreate(((PetscObject)A)->comm,fact);
524: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
525: MatSetType(*fact,((PetscObject)A)->type_name);
526: if (ftype == MAT_FACTOR_LU){
527: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
528: } else {
529: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
530: }
531: (*fact)->factortype = ftype;
532: return(0);
533: }
534: EXTERN_C_END
536: /* ------------------------------------------------------------------*/
539: PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
540: {
541: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
542: PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt;
544: PetscInt m = A->rmap->n,i;
545: PetscBLASInt o = 1,bm = PetscBLASIntCast(m);
548: if (flag & SOR_ZERO_INITIAL_GUESS) {
549: /* this is a hack fix, should have another version without the second BLASdot */
550: VecSet(xx,zero);
551: }
552: VecGetArray(xx,&x);
553: VecGetArray(bb,&b);
554: its = its*lits;
555: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
556: while (its--) {
557: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
558: for (i=0; i<m; i++) {
559: xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
560: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
561: }
562: }
563: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
564: for (i=m-1; i>=0; i--) {
565: xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
566: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
567: }
568: }
569: }
570: VecRestoreArray(bb,&b);
571: VecRestoreArray(xx,&x);
572: return(0);
573: }
575: /* -----------------------------------------------------------------*/
578: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
579: {
580: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
581: PetscScalar *v = mat->v,*x,*y;
583: PetscBLASInt m, n,_One=1;
584: PetscScalar _DOne=1.0,_DZero=0.0;
587: m = PetscBLASIntCast(A->rmap->n);
588: n = PetscBLASIntCast(A->cmap->n);
589: if (!A->rmap->n || !A->cmap->n) return(0);
590: VecGetArray(xx,&x);
591: VecGetArray(yy,&y);
592: BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One);
593: VecRestoreArray(xx,&x);
594: VecRestoreArray(yy,&y);
595: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
596: return(0);
597: }
601: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
602: {
603: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
604: PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
606: PetscBLASInt m, n, _One=1;
609: m = PetscBLASIntCast(A->rmap->n);
610: n = PetscBLASIntCast(A->cmap->n);
611: if (!A->rmap->n || !A->cmap->n) return(0);
612: VecGetArray(xx,&x);
613: VecGetArray(yy,&y);
614: BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One);
615: VecRestoreArray(xx,&x);
616: VecRestoreArray(yy,&y);
617: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
618: return(0);
619: }
623: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
624: {
625: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
626: PetscScalar *v = mat->v,*x,*y,_DOne=1.0;
628: PetscBLASInt m, n, _One=1;
631: m = PetscBLASIntCast(A->rmap->n);
632: n = PetscBLASIntCast(A->cmap->n);
633: if (!A->rmap->n || !A->cmap->n) return(0);
634: if (zz != yy) {VecCopy(zz,yy);}
635: VecGetArray(xx,&x);
636: VecGetArray(yy,&y);
637: BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
638: VecRestoreArray(xx,&x);
639: VecRestoreArray(yy,&y);
640: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
641: return(0);
642: }
646: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
647: {
648: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
649: PetscScalar *v = mat->v,*x,*y;
651: PetscBLASInt m, n, _One=1;
652: PetscScalar _DOne=1.0;
655: m = PetscBLASIntCast(A->rmap->n);
656: n = PetscBLASIntCast(A->cmap->n);
657: if (!A->rmap->n || !A->cmap->n) return(0);
658: if (zz != yy) {VecCopy(zz,yy);}
659: VecGetArray(xx,&x);
660: VecGetArray(yy,&y);
661: BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
662: VecRestoreArray(xx,&x);
663: VecRestoreArray(yy,&y);
664: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
665: return(0);
666: }
668: /* -----------------------------------------------------------------*/
671: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
672: {
673: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
674: PetscScalar *v;
676: PetscInt i;
677:
679: *ncols = A->cmap->n;
680: if (cols) {
681: PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);
682: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
683: }
684: if (vals) {
685: PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);
686: v = mat->v + row;
687: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
688: }
689: return(0);
690: }
694: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
695: {
698: if (cols) {PetscFree(*cols);}
699: if (vals) {PetscFree(*vals); }
700: return(0);
701: }
702: /* ----------------------------------------------------------------*/
705: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
706: {
707: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
708: PetscInt i,j,idx=0;
709:
712: if (!mat->roworiented) {
713: if (addv == INSERT_VALUES) {
714: for (j=0; j<n; j++) {
715: if (indexn[j] < 0) {idx += m; continue;}
716: #if defined(PETSC_USE_DEBUG)
717: 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);
718: #endif
719: for (i=0; i<m; i++) {
720: if (indexm[i] < 0) {idx++; continue;}
721: #if defined(PETSC_USE_DEBUG)
722: 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);
723: #endif
724: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
725: }
726: }
727: } else {
728: for (j=0; j<n; j++) {
729: if (indexn[j] < 0) {idx += m; continue;}
730: #if defined(PETSC_USE_DEBUG)
731: 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);
732: #endif
733: for (i=0; i<m; i++) {
734: if (indexm[i] < 0) {idx++; continue;}
735: #if defined(PETSC_USE_DEBUG)
736: 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);
737: #endif
738: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
739: }
740: }
741: }
742: } else {
743: if (addv == INSERT_VALUES) {
744: for (i=0; i<m; i++) {
745: if (indexm[i] < 0) { idx += n; continue;}
746: #if defined(PETSC_USE_DEBUG)
747: 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);
748: #endif
749: for (j=0; j<n; j++) {
750: if (indexn[j] < 0) { idx++; continue;}
751: #if defined(PETSC_USE_DEBUG)
752: 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);
753: #endif
754: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
755: }
756: }
757: } else {
758: for (i=0; i<m; i++) {
759: if (indexm[i] < 0) { idx += n; continue;}
760: #if defined(PETSC_USE_DEBUG)
761: 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);
762: #endif
763: for (j=0; j<n; j++) {
764: if (indexn[j] < 0) { idx++; continue;}
765: #if defined(PETSC_USE_DEBUG)
766: 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);
767: #endif
768: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
769: }
770: }
771: }
772: }
773: return(0);
774: }
778: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
779: {
780: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
781: PetscInt i,j;
784: /* row-oriented output */
785: for (i=0; i<m; i++) {
786: if (indexm[i] < 0) {v += n;continue;}
787: 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);
788: for (j=0; j<n; j++) {
789: if (indexn[j] < 0) {v++; continue;}
790: 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);
791: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
792: }
793: }
794: return(0);
795: }
797: /* -----------------------------------------------------------------*/
801: PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
802: {
803: Mat_SeqDense *a;
805: PetscInt *scols,i,j,nz,header[4];
806: int fd;
807: PetscMPIInt size;
808: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
809: PetscScalar *vals,*svals,*v,*w;
810: MPI_Comm comm = ((PetscObject)viewer)->comm;
813: MPI_Comm_size(comm,&size);
814: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
815: PetscViewerBinaryGetDescriptor(viewer,&fd);
816: PetscBinaryRead(fd,header,4,PETSC_INT);
817: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
818: M = header[1]; N = header[2]; nz = header[3];
820: /* set global size if not set already*/
821: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
822: MatSetSizes(newmat,M,N,M,N);
823: } else {
824: /* if sizes and type are already set, check if the vector global sizes are correct */
825: MatGetSize(newmat,&grows,&gcols);
826: 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);
827: }
828: MatSeqDenseSetPreallocation(newmat,PETSC_NULL);
829:
830: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
831: a = (Mat_SeqDense*)newmat->data;
832: v = a->v;
833: /* Allocate some temp space to read in the values and then flip them
834: from row major to column major */
835: PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);
836: /* read in nonzero values */
837: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
838: /* now flip the values and store them in the matrix*/
839: for (j=0; j<N; j++) {
840: for (i=0; i<M; i++) {
841: *v++ =w[i*N+j];
842: }
843: }
844: PetscFree(w);
845: } else {
846: /* read row lengths */
847: PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);
848: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
850: a = (Mat_SeqDense*)newmat->data;
851: v = a->v;
853: /* read column indices and nonzeros */
854: PetscMalloc((nz+1)*sizeof(PetscInt),&scols);
855: cols = scols;
856: PetscBinaryRead(fd,cols,nz,PETSC_INT);
857: PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);
858: vals = svals;
859: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
861: /* insert into matrix */
862: for (i=0; i<M; i++) {
863: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
864: svals += rowlengths[i]; scols += rowlengths[i];
865: }
866: PetscFree(vals);
867: PetscFree(cols);
868: PetscFree(rowlengths);
869: }
870: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
871: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
873: return(0);
874: }
878: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
879: {
880: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
881: PetscErrorCode ierr;
882: PetscInt i,j;
883: const char *name;
884: PetscScalar *v;
885: PetscViewerFormat format;
886: #if defined(PETSC_USE_COMPLEX)
887: PetscBool allreal = PETSC_TRUE;
888: #endif
891: PetscViewerGetFormat(viewer,&format);
892: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
893: return(0); /* do nothing for now */
894: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
895: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
896: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
897: for (i=0; i<A->rmap->n; i++) {
898: v = a->v + i;
899: PetscViewerASCIIPrintf(viewer,"row %D:",i);
900: for (j=0; j<A->cmap->n; j++) {
901: #if defined(PETSC_USE_COMPLEX)
902: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
903: PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));
904: } else if (PetscRealPart(*v)) {
905: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));
906: }
907: #else
908: if (*v) {
909: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);
910: }
911: #endif
912: v += a->lda;
913: }
914: PetscViewerASCIIPrintf(viewer,"\n");
915: }
916: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
917: } else {
918: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
919: #if defined(PETSC_USE_COMPLEX)
920: /* determine if matrix has all real values */
921: v = a->v;
922: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
923: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;}
924: }
925: #endif
926: if (format == PETSC_VIEWER_ASCII_MATLAB) {
927: PetscObjectGetName((PetscObject)A,&name);
928: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
929: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
930: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
931: } else {
932: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
933: }
935: for (i=0; i<A->rmap->n; i++) {
936: v = a->v + i;
937: for (j=0; j<A->cmap->n; j++) {
938: #if defined(PETSC_USE_COMPLEX)
939: if (allreal) {
940: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
941: } else {
942: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
943: }
944: #else
945: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
946: #endif
947: v += a->lda;
948: }
949: PetscViewerASCIIPrintf(viewer,"\n");
950: }
951: if (format == PETSC_VIEWER_ASCII_MATLAB) {
952: PetscViewerASCIIPrintf(viewer,"];\n");
953: }
954: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
955: }
956: PetscViewerFlush(viewer);
957: return(0);
958: }
962: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
963: {
964: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
965: PetscErrorCode ierr;
966: int fd;
967: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
968: PetscScalar *v,*anonz,*vals;
969: PetscViewerFormat format;
970:
972: PetscViewerBinaryGetDescriptor(viewer,&fd);
974: PetscViewerGetFormat(viewer,&format);
975: if (format == PETSC_VIEWER_NATIVE) {
976: /* store the matrix as a dense matrix */
977: PetscMalloc(4*sizeof(PetscInt),&col_lens);
978: col_lens[0] = MAT_FILE_CLASSID;
979: col_lens[1] = m;
980: col_lens[2] = n;
981: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
982: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
983: PetscFree(col_lens);
985: /* write out matrix, by rows */
986: PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);
987: v = a->v;
988: for (j=0; j<n; j++) {
989: for (i=0; i<m; i++) {
990: vals[j + i*n] = *v++;
991: }
992: }
993: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
994: PetscFree(vals);
995: } else {
996: PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);
997: col_lens[0] = MAT_FILE_CLASSID;
998: col_lens[1] = m;
999: col_lens[2] = n;
1000: col_lens[3] = nz;
1002: /* store lengths of each row and write (including header) to file */
1003: for (i=0; i<m; i++) col_lens[4+i] = n;
1004: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1006: /* Possibly should write in smaller increments, not whole matrix at once? */
1007: /* store column indices (zero start index) */
1008: ict = 0;
1009: for (i=0; i<m; i++) {
1010: for (j=0; j<n; j++) col_lens[ict++] = j;
1011: }
1012: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1013: PetscFree(col_lens);
1015: /* store nonzero values */
1016: PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);
1017: ict = 0;
1018: for (i=0; i<m; i++) {
1019: v = a->v + i;
1020: for (j=0; j<n; j++) {
1021: anonz[ict++] = *v; v += a->lda;
1022: }
1023: }
1024: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1025: PetscFree(anonz);
1026: }
1027: return(0);
1028: }
1032: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1033: {
1034: Mat A = (Mat) Aa;
1035: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1036: PetscErrorCode ierr;
1037: PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j;
1038: PetscScalar *v = a->v;
1039: PetscViewer viewer;
1040: PetscDraw popup;
1041: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
1042: PetscViewerFormat format;
1046: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1047: PetscViewerGetFormat(viewer,&format);
1048: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1050: /* Loop over matrix elements drawing boxes */
1051: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1052: /* Blue for negative and Red for positive */
1053: color = PETSC_DRAW_BLUE;
1054: for(j = 0; j < n; j++) {
1055: x_l = j;
1056: x_r = x_l + 1.0;
1057: for(i = 0; i < m; i++) {
1058: y_l = m - i - 1.0;
1059: y_r = y_l + 1.0;
1060: #if defined(PETSC_USE_COMPLEX)
1061: if (PetscRealPart(v[j*m+i]) > 0.) {
1062: color = PETSC_DRAW_RED;
1063: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1064: color = PETSC_DRAW_BLUE;
1065: } else {
1066: continue;
1067: }
1068: #else
1069: if (v[j*m+i] > 0.) {
1070: color = PETSC_DRAW_RED;
1071: } else if (v[j*m+i] < 0.) {
1072: color = PETSC_DRAW_BLUE;
1073: } else {
1074: continue;
1075: }
1076: #endif
1077: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1078: }
1079: }
1080: } else {
1081: /* use contour shading to indicate magnitude of values */
1082: /* first determine max of all nonzero values */
1083: for(i = 0; i < m*n; i++) {
1084: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1085: }
1086: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
1087: PetscDrawGetPopup(draw,&popup);
1088: if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
1089: for(j = 0; j < n; j++) {
1090: x_l = j;
1091: x_r = x_l + 1.0;
1092: for(i = 0; i < m; i++) {
1093: y_l = m - i - 1.0;
1094: y_r = y_l + 1.0;
1095: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
1096: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1097: }
1098: }
1099: }
1100: return(0);
1101: }
1105: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1106: {
1107: PetscDraw draw;
1108: PetscBool isnull;
1109: PetscReal xr,yr,xl,yl,h,w;
1113: PetscViewerDrawGetDraw(viewer,0,&draw);
1114: PetscDrawIsNull(draw,&isnull);
1115: if (isnull) return(0);
1117: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1118: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1119: xr += w; yr += h; xl = -w; yl = -h;
1120: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1121: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1122: PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
1123: return(0);
1124: }
1128: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1129: {
1131: PetscBool iascii,isbinary,isdraw;
1134: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1135: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1136: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1138: if (iascii) {
1139: MatView_SeqDense_ASCII(A,viewer);
1140: } else if (isbinary) {
1141: MatView_SeqDense_Binary(A,viewer);
1142: } else if (isdraw) {
1143: MatView_SeqDense_Draw(A,viewer);
1144: } else {
1145: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name);
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: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);
1167: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);
1168: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);
1169: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);
1170: return(0);
1171: }
1175: PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1176: {
1177: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1179: PetscInt k,j,m,n,M;
1180: PetscScalar *v,tmp;
1183: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1184: if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */
1185: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1186: else {
1187: for (j=0; j<m; j++) {
1188: for (k=0; k<j; k++) {
1189: tmp = v[j + k*M];
1190: v[j + k*M] = v[k + j*M];
1191: v[k + j*M] = tmp;
1192: }
1193: }
1194: }
1195: } else { /* out-of-place transpose */
1196: Mat tmat;
1197: Mat_SeqDense *tmatd;
1198: PetscScalar *v2;
1200: if (reuse == MAT_INITIAL_MATRIX) {
1201: MatCreate(((PetscObject)A)->comm,&tmat);
1202: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1203: MatSetType(tmat,((PetscObject)A)->type_name);
1204: MatSeqDenseSetPreallocation(tmat,PETSC_NULL);
1205: } else {
1206: tmat = *matout;
1207: }
1208: tmatd = (Mat_SeqDense*)tmat->data;
1209: v = mat->v; v2 = tmatd->v;
1210: for (j=0; j<n; j++) {
1211: for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
1212: }
1213: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1214: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1215: *matout = tmat;
1216: }
1217: return(0);
1218: }
1222: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1223: {
1224: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1225: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1226: PetscInt i,j;
1227: PetscScalar *v1 = mat1->v,*v2 = mat2->v;
1230: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1231: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1232: for (i=0; i<A1->rmap->n; i++) {
1233: v1 = mat1->v+i; v2 = mat2->v+i;
1234: for (j=0; j<A1->cmap->n; j++) {
1235: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1236: v1 += mat1->lda; v2 += mat2->lda;
1237: }
1238: }
1239: *flg = PETSC_TRUE;
1240: return(0);
1241: }
1245: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1246: {
1247: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1249: PetscInt i,n,len;
1250: PetscScalar *x,zero = 0.0;
1253: VecSet(v,zero);
1254: VecGetSize(v,&n);
1255: VecGetArray(v,&x);
1256: len = PetscMin(A->rmap->n,A->cmap->n);
1257: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1258: for (i=0; i<len; i++) {
1259: x[i] = mat->v[i*mat->lda + i];
1260: }
1261: VecRestoreArray(v,&x);
1262: return(0);
1263: }
1267: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1268: {
1269: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1270: PetscScalar *l,*r,x,*v;
1272: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1275: if (ll) {
1276: VecGetSize(ll,&m);
1277: VecGetArray(ll,&l);
1278: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1279: for (i=0; i<m; i++) {
1280: x = l[i];
1281: v = mat->v + i;
1282: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1283: }
1284: VecRestoreArray(ll,&l);
1285: PetscLogFlops(n*m);
1286: }
1287: if (rr) {
1288: VecGetSize(rr,&n);
1289: VecGetArray(rr,&r);
1290: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1291: for (i=0; i<n; i++) {
1292: x = r[i];
1293: v = mat->v + i*m;
1294: for (j=0; j<m; j++) { (*v++) *= x;}
1295: }
1296: VecRestoreArray(rr,&r);
1297: PetscLogFlops(n*m);
1298: }
1299: return(0);
1300: }
1304: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1305: {
1306: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1307: PetscScalar *v = mat->v;
1308: PetscReal sum = 0.0;
1309: PetscInt lda=mat->lda,m=A->rmap->n,i,j;
1313: if (type == NORM_FROBENIUS) {
1314: if (lda>m) {
1315: for (j=0; j<A->cmap->n; j++) {
1316: v = mat->v+j*lda;
1317: for (i=0; i<m; i++) {
1318: #if defined(PETSC_USE_COMPLEX)
1319: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1320: #else
1321: sum += (*v)*(*v); v++;
1322: #endif
1323: }
1324: }
1325: } else {
1326: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1327: #if defined(PETSC_USE_COMPLEX)
1328: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1329: #else
1330: sum += (*v)*(*v); v++;
1331: #endif
1332: }
1333: }
1334: *nrm = PetscSqrtReal(sum);
1335: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1336: } else if (type == NORM_1) {
1337: *nrm = 0.0;
1338: for (j=0; j<A->cmap->n; j++) {
1339: v = mat->v + j*mat->lda;
1340: sum = 0.0;
1341: for (i=0; i<A->rmap->n; i++) {
1342: sum += PetscAbsScalar(*v); v++;
1343: }
1344: if (sum > *nrm) *nrm = sum;
1345: }
1346: PetscLogFlops(A->cmap->n*A->rmap->n);
1347: } else if (type == NORM_INFINITY) {
1348: *nrm = 0.0;
1349: for (j=0; j<A->rmap->n; j++) {
1350: v = mat->v + j;
1351: sum = 0.0;
1352: for (i=0; i<A->cmap->n; i++) {
1353: sum += PetscAbsScalar(*v); v += mat->lda;
1354: }
1355: if (sum > *nrm) *nrm = sum;
1356: }
1357: PetscLogFlops(A->cmap->n*A->rmap->n);
1358: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1359: return(0);
1360: }
1364: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1365: {
1366: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1368:
1370: switch (op) {
1371: case MAT_ROW_ORIENTED:
1372: aij->roworiented = flg;
1373: break;
1374: case MAT_NEW_NONZERO_LOCATIONS:
1375: case MAT_NEW_NONZERO_LOCATION_ERR:
1376: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1377: case MAT_NEW_DIAGONALS:
1378: case MAT_KEEP_NONZERO_PATTERN:
1379: case MAT_IGNORE_OFF_PROC_ENTRIES:
1380: case MAT_USE_HASH_TABLE:
1381: case MAT_SYMMETRIC:
1382: case MAT_STRUCTURALLY_SYMMETRIC:
1383: case MAT_HERMITIAN:
1384: case MAT_SYMMETRY_ETERNAL:
1385: case MAT_IGNORE_LOWER_TRIANGULAR:
1386: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1387: break;
1388: default:
1389: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1390: }
1391: return(0);
1392: }
1396: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1397: {
1398: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1400: PetscInt lda=l->lda,m=A->rmap->n,j;
1403: if (lda>m) {
1404: for (j=0; j<A->cmap->n; j++) {
1405: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1406: }
1407: } else {
1408: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1409: }
1410: return(0);
1411: }
1415: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1416: {
1417: PetscErrorCode ierr;
1418: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1419: PetscInt m = l->lda, n = A->cmap->n, i,j;
1420: PetscScalar *slot,*bb;
1421: const PetscScalar *xx;
1424: #if defined(PETSC_USE_DEBUG)
1425: for (i=0; i<N; i++) {
1426: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1427: 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);
1428: }
1429: #endif
1431: /* fix right hand side if needed */
1432: if (x && b) {
1433: VecGetArrayRead(x,&xx);
1434: VecGetArray(b,&bb);
1435: for (i=0; i<N; i++) {
1436: bb[rows[i]] = diag*xx[rows[i]];
1437: }
1438: VecRestoreArrayRead(x,&xx);
1439: VecRestoreArray(b,&bb);
1440: }
1442: for (i=0; i<N; i++) {
1443: slot = l->v + rows[i];
1444: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1445: }
1446: if (diag != 0.0) {
1447: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1448: for (i=0; i<N; i++) {
1449: slot = l->v + (m+1)*rows[i];
1450: *slot = diag;
1451: }
1452: }
1453: return(0);
1454: }
1458: PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[])
1459: {
1460: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1463: 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");
1464: *array = mat->v;
1465: return(0);
1466: }
1470: PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1471: {
1473: *array = 0; /* user cannot accidently use the array later */
1474: return(0);
1475: }
1479: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1480: {
1481: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1483: PetscInt i,j,nrows,ncols;
1484: const PetscInt *irow,*icol;
1485: PetscScalar *av,*bv,*v = mat->v;
1486: Mat newmat;
1489: ISGetIndices(isrow,&irow);
1490: ISGetIndices(iscol,&icol);
1491: ISGetLocalSize(isrow,&nrows);
1492: ISGetLocalSize(iscol,&ncols);
1493:
1494: /* Check submatrixcall */
1495: if (scall == MAT_REUSE_MATRIX) {
1496: PetscInt n_cols,n_rows;
1497: MatGetSize(*B,&n_rows,&n_cols);
1498: if (n_rows != nrows || n_cols != ncols) {
1499: /* resize the result matrix to match number of requested rows/columns */
1500: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1501: }
1502: newmat = *B;
1503: } else {
1504: /* Create and fill new matrix */
1505: MatCreate(((PetscObject)A)->comm,&newmat);
1506: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1507: MatSetType(newmat,((PetscObject)A)->type_name);
1508: MatSeqDenseSetPreallocation(newmat,PETSC_NULL);
1509: }
1511: /* Now extract the data pointers and do the copy,column at a time */
1512: bv = ((Mat_SeqDense*)newmat->data)->v;
1513:
1514: for (i=0; i<ncols; i++) {
1515: av = v + mat->lda*icol[i];
1516: for (j=0; j<nrows; j++) {
1517: *bv++ = av[irow[j]];
1518: }
1519: }
1521: /* Assemble the matrices so that the correct flags are set */
1522: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1523: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1525: /* Free work space */
1526: ISRestoreIndices(isrow,&irow);
1527: ISRestoreIndices(iscol,&icol);
1528: *B = newmat;
1529: return(0);
1530: }
1534: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1535: {
1537: PetscInt i;
1540: if (scall == MAT_INITIAL_MATRIX) {
1541: PetscMalloc((n+1)*sizeof(Mat),B);
1542: }
1544: for (i=0; i<n; i++) {
1545: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1546: }
1547: return(0);
1548: }
1552: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1553: {
1555: return(0);
1556: }
1560: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1561: {
1563: return(0);
1564: }
1568: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1569: {
1570: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data;
1572: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1575: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1576: if (A->ops->copy != B->ops->copy) {
1577: MatCopy_Basic(A,B,str);
1578: return(0);
1579: }
1580: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1581: if (lda1>m || lda2>m) {
1582: for (j=0; j<n; j++) {
1583: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1584: }
1585: } else {
1586: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1587: }
1588: return(0);
1589: }
1593: PetscErrorCode MatSetUp_SeqDense(Mat A)
1594: {
1598: MatSeqDenseSetPreallocation(A,0);
1599: return(0);
1600: }
1604: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1605: {
1606: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1607: PetscInt i,nz = A->rmap->n*A->cmap->n;
1608: PetscScalar *aa = a->v;
1611: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1612: return(0);
1613: }
1617: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1618: {
1619: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1620: PetscInt i,nz = A->rmap->n*A->cmap->n;
1621: PetscScalar *aa = a->v;
1624: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1625: return(0);
1626: }
1630: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1631: {
1632: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1633: PetscInt i,nz = A->rmap->n*A->cmap->n;
1634: PetscScalar *aa = a->v;
1637: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1638: return(0);
1639: }
1641: /* ----------------------------------------------------------------*/
1644: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1645: {
1649: if (scall == MAT_INITIAL_MATRIX){
1650: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1651: }
1652: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1653: return(0);
1654: }
1658: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1659: {
1661: PetscInt m=A->rmap->n,n=B->cmap->n;
1662: Mat Cmat;
1665: 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);
1666: MatCreate(PETSC_COMM_SELF,&Cmat);
1667: MatSetSizes(Cmat,m,n,m,n);
1668: MatSetType(Cmat,MATSEQDENSE);
1669: MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);
1670: Cmat->assembled = PETSC_TRUE;
1671: Cmat->ops->matmult = MatMatMult_SeqDense_SeqDense;
1672: *C = Cmat;
1673: return(0);
1674: }
1678: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1679: {
1680: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1681: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1682: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1683: PetscBLASInt m,n,k;
1684: PetscScalar _DOne=1.0,_DZero=0.0;
1687: m = PetscBLASIntCast(A->rmap->n);
1688: n = PetscBLASIntCast(B->cmap->n);
1689: k = PetscBLASIntCast(A->cmap->n);
1690: BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda);
1691: return(0);
1692: }
1696: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1697: {
1701: if (scall == MAT_INITIAL_MATRIX){
1702: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1703: }
1704: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1705: return(0);
1706: }
1710: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1711: {
1713: PetscInt m=A->cmap->n,n=B->cmap->n;
1714: Mat Cmat;
1717: 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);
1718: MatCreate(PETSC_COMM_SELF,&Cmat);
1719: MatSetSizes(Cmat,m,n,m,n);
1720: MatSetType(Cmat,MATSEQDENSE);
1721: MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);
1722: Cmat->assembled = PETSC_TRUE;
1723: *C = Cmat;
1724: return(0);
1725: }
1729: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1730: {
1731: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1732: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1733: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1734: PetscBLASInt m,n,k;
1735: PetscScalar _DOne=1.0,_DZero=0.0;
1738: m = PetscBLASIntCast(A->cmap->n);
1739: n = PetscBLASIntCast(B->cmap->n);
1740: k = PetscBLASIntCast(A->rmap->n);
1741: /*
1742: Note the m and n arguments below are the number rows and columns of A', not A!
1743: */
1744: BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda);
1745: return(0);
1746: }
1750: PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
1751: {
1752: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1754: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1755: PetscScalar *x;
1756: MatScalar *aa = a->v;
1759: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1761: VecSet(v,0.0);
1762: VecGetArray(v,&x);
1763: VecGetLocalSize(v,&p);
1764: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1765: for (i=0; i<m; i++) {
1766: x[i] = aa[i]; if (idx) idx[i] = 0;
1767: for (j=1; j<n; j++){
1768: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1769: }
1770: }
1771: VecRestoreArray(v,&x);
1772: return(0);
1773: }
1777: PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
1778: {
1779: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1781: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1782: PetscScalar *x;
1783: PetscReal atmp;
1784: MatScalar *aa = a->v;
1787: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1789: VecSet(v,0.0);
1790: VecGetArray(v,&x);
1791: VecGetLocalSize(v,&p);
1792: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1793: for (i=0; i<m; i++) {
1794: x[i] = PetscAbsScalar(aa[i]);
1795: for (j=1; j<n; j++){
1796: atmp = PetscAbsScalar(aa[i+m*j]);
1797: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
1798: }
1799: }
1800: VecRestoreArray(v,&x);
1801: return(0);
1802: }
1806: PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
1807: {
1808: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1810: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
1811: PetscScalar *x;
1812: MatScalar *aa = a->v;
1815: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1817: VecSet(v,0.0);
1818: VecGetArray(v,&x);
1819: VecGetLocalSize(v,&p);
1820: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1821: for (i=0; i<m; i++) {
1822: x[i] = aa[i]; if (idx) idx[i] = 0;
1823: for (j=1; j<n; j++){
1824: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
1825: }
1826: }
1827: VecRestoreArray(v,&x);
1828: return(0);
1829: }
1833: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
1834: {
1835: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1837: PetscScalar *x;
1840: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1842: VecGetArray(v,&x);
1843: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
1844: VecRestoreArray(v,&x);
1845: return(0);
1846: }
1851: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
1852: {
1854: PetscInt i,j,m,n;
1855: PetscScalar *a;
1858: MatGetSize(A,&m,&n);
1859: PetscMemzero(norms,n*sizeof(PetscReal));
1860: MatGetArray(A,&a);
1861: if (type == NORM_2) {
1862: for (i=0; i<n; i++ ){
1863: for (j=0; j<m; j++) {
1864: norms[i] += PetscAbsScalar(a[j]*a[j]);
1865: }
1866: a += m;
1867: }
1868: } else if (type == NORM_1) {
1869: for (i=0; i<n; i++ ){
1870: for (j=0; j<m; j++) {
1871: norms[i] += PetscAbsScalar(a[j]);
1872: }
1873: a += m;
1874: }
1875: } else if (type == NORM_INFINITY) {
1876: for (i=0; i<n; i++ ){
1877: for (j=0; j<m; j++) {
1878: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
1879: }
1880: a += m;
1881: }
1882: } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType");
1883: if (type == NORM_2) {
1884: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
1885: }
1886: return(0);
1887: }
1889: /* -------------------------------------------------------------------*/
1890: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
1891: MatGetRow_SeqDense,
1892: MatRestoreRow_SeqDense,
1893: MatMult_SeqDense,
1894: /* 4*/ MatMultAdd_SeqDense,
1895: MatMultTranspose_SeqDense,
1896: MatMultTransposeAdd_SeqDense,
1897: 0,
1898: 0,
1899: 0,
1900: /*10*/ 0,
1901: MatLUFactor_SeqDense,
1902: MatCholeskyFactor_SeqDense,
1903: MatSOR_SeqDense,
1904: MatTranspose_SeqDense,
1905: /*15*/ MatGetInfo_SeqDense,
1906: MatEqual_SeqDense,
1907: MatGetDiagonal_SeqDense,
1908: MatDiagonalScale_SeqDense,
1909: MatNorm_SeqDense,
1910: /*20*/ MatAssemblyBegin_SeqDense,
1911: MatAssemblyEnd_SeqDense,
1912: MatSetOption_SeqDense,
1913: MatZeroEntries_SeqDense,
1914: /*24*/ MatZeroRows_SeqDense,
1915: 0,
1916: 0,
1917: 0,
1918: 0,
1919: /*29*/ MatSetUp_SeqDense,
1920: 0,
1921: 0,
1922: MatGetArray_SeqDense,
1923: MatRestoreArray_SeqDense,
1924: /*34*/ MatDuplicate_SeqDense,
1925: 0,
1926: 0,
1927: 0,
1928: 0,
1929: /*39*/ MatAXPY_SeqDense,
1930: MatGetSubMatrices_SeqDense,
1931: 0,
1932: MatGetValues_SeqDense,
1933: MatCopy_SeqDense,
1934: /*44*/ MatGetRowMax_SeqDense,
1935: MatScale_SeqDense,
1936: 0,
1937: 0,
1938: 0,
1939: /*49*/ 0,
1940: 0,
1941: 0,
1942: 0,
1943: 0,
1944: /*54*/ 0,
1945: 0,
1946: 0,
1947: 0,
1948: 0,
1949: /*59*/ 0,
1950: MatDestroy_SeqDense,
1951: MatView_SeqDense,
1952: 0,
1953: 0,
1954: /*64*/ 0,
1955: 0,
1956: 0,
1957: 0,
1958: 0,
1959: /*69*/ MatGetRowMaxAbs_SeqDense,
1960: 0,
1961: 0,
1962: 0,
1963: 0,
1964: /*74*/ 0,
1965: 0,
1966: 0,
1967: 0,
1968: 0,
1969: /*79*/ 0,
1970: 0,
1971: 0,
1972: 0,
1973: /*83*/ MatLoad_SeqDense,
1974: 0,
1975: MatIsHermitian_SeqDense,
1976: 0,
1977: 0,
1978: 0,
1979: /*89*/ MatMatMult_SeqDense_SeqDense,
1980: MatMatMultSymbolic_SeqDense_SeqDense,
1981: MatMatMultNumeric_SeqDense_SeqDense,
1982: 0,
1983: 0,
1984: /*94*/ 0,
1985: 0,
1986: 0,
1987: 0,
1988: 0,
1989: /*99*/ 0,
1990: 0,
1991: 0,
1992: MatConjugate_SeqDense,
1993: 0,
1994: /*104*/0,
1995: MatRealPart_SeqDense,
1996: MatImaginaryPart_SeqDense,
1997: 0,
1998: 0,
1999: /*109*/MatMatSolve_SeqDense,
2000: 0,
2001: MatGetRowMin_SeqDense,
2002: MatGetColumnVector_SeqDense,
2003: 0,
2004: /*114*/0,
2005: 0,
2006: 0,
2007: 0,
2008: 0,
2009: /*119*/0,
2010: 0,
2011: 0,
2012: 0,
2013: 0,
2014: /*124*/0,
2015: MatGetColumnNorms_SeqDense,
2016: 0,
2017: 0,
2018: 0,
2019: /*129*/0,
2020: MatTransposeMatMult_SeqDense_SeqDense,
2021: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2022: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2023: };
2027: /*@C
2028: MatCreateSeqDense - Creates a sequential dense matrix that
2029: is stored in column major order (the usual Fortran 77 manner). Many
2030: of the matrix operations use the BLAS and LAPACK routines.
2032: Collective on MPI_Comm
2034: Input Parameters:
2035: + comm - MPI communicator, set to PETSC_COMM_SELF
2036: . m - number of rows
2037: . n - number of columns
2038: - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc
2039: to control all matrix memory allocation.
2041: Output Parameter:
2042: . A - the matrix
2044: Notes:
2045: The data input variable is intended primarily for Fortran programmers
2046: who wish to allocate their own matrix memory space. Most users should
2047: set data=PETSC_NULL.
2049: Level: intermediate
2051: .keywords: dense, matrix, LAPACK, BLAS
2053: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2054: @*/
2055: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2056: {
2060: MatCreate(comm,A);
2061: MatSetSizes(*A,m,n,m,n);
2062: MatSetType(*A,MATSEQDENSE);
2063: MatSeqDenseSetPreallocation(*A,data);
2064: return(0);
2065: }
2069: /*@C
2070: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2072: Collective on MPI_Comm
2074: Input Parameters:
2075: + A - the matrix
2076: - data - the array (or PETSC_NULL)
2078: Notes:
2079: The data input variable is intended primarily for Fortran programmers
2080: who wish to allocate their own matrix memory space. Most users should
2081: need not call this routine.
2083: Level: intermediate
2085: .keywords: dense, matrix, LAPACK, BLAS
2087: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2089: @*/
2090: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2091: {
2095: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2096: return(0);
2097: }
2099: EXTERN_C_BEGIN
2102: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2103: {
2104: Mat_SeqDense *b;
2108: B->preallocated = PETSC_TRUE;
2110: PetscLayoutSetUp(B->rmap);
2111: PetscLayoutSetUp(B->cmap);
2113: b = (Mat_SeqDense*)B->data;
2114: b->Mmax = B->rmap->n;
2115: b->Nmax = B->cmap->n;
2116: if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2118: if (!data) { /* petsc-allocated storage */
2119: if (!b->user_alloc) { PetscFree(b->v); }
2120: PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);
2121: PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));
2122: PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));
2123: b->user_alloc = PETSC_FALSE;
2124: } else { /* user-allocated storage */
2125: if (!b->user_alloc) { PetscFree(b->v); }
2126: b->v = data;
2127: b->user_alloc = PETSC_TRUE;
2128: }
2129: B->assembled = PETSC_TRUE;
2130: return(0);
2131: }
2132: EXTERN_C_END
2136: /*@C
2137: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2139: Input parameter:
2140: + A - the matrix
2141: - lda - the leading dimension
2143: Notes:
2144: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2145: it asserts that the preallocation has a leading dimension (the LDA parameter
2146: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2148: Level: intermediate
2150: .keywords: dense, matrix, LAPACK, BLAS
2152: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2154: @*/
2155: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2156: {
2157: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2160: 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);
2161: b->lda = lda;
2162: b->changelda = PETSC_FALSE;
2163: b->Mmax = PetscMax(b->Mmax,lda);
2164: return(0);
2165: }
2167: /*MC
2168: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2170: Options Database Keys:
2171: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2173: Level: beginner
2175: .seealso: MatCreateSeqDense()
2177: M*/
2179: EXTERN_C_BEGIN
2182: PetscErrorCode MatCreate_SeqDense(Mat B)
2183: {
2184: Mat_SeqDense *b;
2186: PetscMPIInt size;
2189: MPI_Comm_size(((PetscObject)B)->comm,&size);
2190: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2192: PetscNewLog(B,Mat_SeqDense,&b);
2193: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2194: B->data = (void*)b;
2196: b->pivots = 0;
2197: b->roworiented = PETSC_TRUE;
2198: b->v = 0;
2199: b->changelda = PETSC_FALSE;
2202: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_seqaij_C","MatConvert_SeqDense_SeqAIJ",MatConvert_SeqDense_SeqAIJ);
2203: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C",
2204: "MatGetFactor_seqdense_petsc",
2205: MatGetFactor_seqdense_petsc);
2206: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C",
2207: "MatSeqDenseSetPreallocation_SeqDense",
2208: MatSeqDenseSetPreallocation_SeqDense);
2210: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C",
2211: "MatMatMult_SeqAIJ_SeqDense",
2212: MatMatMult_SeqAIJ_SeqDense);
2215: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",
2216: "MatMatMultSymbolic_SeqAIJ_SeqDense",
2217: MatMatMultSymbolic_SeqAIJ_SeqDense);
2218: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",
2219: "MatMatMultNumeric_SeqAIJ_SeqDense",
2220: MatMatMultNumeric_SeqAIJ_SeqDense);
2221: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2222: return(0);
2223: }
2224: EXTERN_C_END