Actual source code: dense.c
1: #define PETSCMAT_DLL
3: /*
4: Defines the basic matrix operations for sequential dense.
5: */
7: #include src/mat/impls/dense/seq/dense.h
8: #include petscblaslapack.h
12: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
13: {
14: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
15: PetscScalar oalpha = alpha;
16: PetscInt j;
17: PetscBLASInt N = (PetscBLASInt)X->rmap.n*X->cmap.n,m=(PetscBLASInt)X->rmap.n,ldax = x->lda,lday=y->lda,one = 1;
21: if (ldax>m || lday>m) {
22: for (j=0; j<X->cmap.n; j++) {
23: BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one);
24: }
25: } else {
26: BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one);
27: }
28: PetscLogFlops(2*N-1);
29: return(0);
30: }
34: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
35: {
36: PetscInt N = A->rmap.n*A->cmap.n;
39: info->rows_global = (double)A->rmap.n;
40: info->columns_global = (double)A->cmap.n;
41: info->rows_local = (double)A->rmap.n;
42: info->columns_local = (double)A->cmap.n;
43: info->block_size = 1.0;
44: info->nz_allocated = (double)N;
45: info->nz_used = (double)N;
46: info->nz_unneeded = (double)0;
47: info->assemblies = (double)A->num_ass;
48: info->mallocs = 0;
49: info->memory = A->mem;
50: info->fill_ratio_given = 0;
51: info->fill_ratio_needed = 0;
52: info->factor_mallocs = 0;
53: return(0);
54: }
58: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
59: {
60: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
61: PetscScalar oalpha = alpha;
62: PetscBLASInt one = 1,lda = a->lda,j,nz;
66: if (lda>A->rmap.n) {
67: nz = (PetscBLASInt)A->rmap.n;
68: for (j=0; j<A->cmap.n; j++) {
69: BLASscal_(&nz,&oalpha,a->v+j*lda,&one);
70: }
71: } else {
72: nz = (PetscBLASInt)A->rmap.n*A->cmap.n;
73: BLASscal_(&nz,&oalpha,a->v,&one);
74: }
75: PetscLogFlops(nz);
76: return(0);
77: }
78:
79: /* ---------------------------------------------------------------*/
80: /* COMMENT: I have chosen to hide row permutation in the pivots,
81: rather than put it in the Mat->row slot.*/
84: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo)
85: {
86: #if defined(PETSC_MISSING_LAPACK_GETRF)
88: SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
89: #else
90: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
92: PetscBLASInt n = (PetscBLASInt)A->cmap.n,m = (PetscBLASInt)A->rmap.n,info;
95: if (!mat->pivots) {
96: PetscMalloc((A->rmap.n+1)*sizeof(PetscBLASInt),&mat->pivots);
97: PetscLogObjectMemory(A,A->rmap.n*sizeof(PetscBLASInt));
98: }
99: A->factor = FACTOR_LU;
100: if (!A->rmap.n || !A->cmap.n) return(0);
101: LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info);
102: if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization");
103: if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
104: PetscLogFlops((2*A->cmap.n*A->cmap.n*A->cmap.n)/3);
105: #endif
106: return(0);
107: }
111: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
112: {
113: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
115: PetscInt lda = (PetscInt)mat->lda,j,m;
116: Mat newi;
119: MatCreate(A->comm,&newi);
120: MatSetSizes(newi,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
121: MatSetType(newi,A->type_name);
122: MatSeqDenseSetPreallocation(newi,PETSC_NULL);
123: if (cpvalues == MAT_COPY_VALUES) {
124: l = (Mat_SeqDense*)newi->data;
125: if (lda>A->rmap.n) {
126: m = A->rmap.n;
127: for (j=0; j<A->cmap.n; j++) {
128: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
129: }
130: } else {
131: PetscMemcpy(l->v,mat->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));
132: }
133: }
134: newi->assembled = PETSC_TRUE;
135: *newmat = newi;
136: return(0);
137: }
141: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatFactorInfo *info,Mat *fact)
142: {
146: MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);
147: return(0);
148: }
152: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact)
153: {
154: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l = (Mat_SeqDense*)(*fact)->data;
156: PetscInt lda1=mat->lda,lda2=l->lda, m=A->rmap.n,n=A->cmap.n, j;
157: MatFactorInfo info;
160: /* copy the numerical values */
161: if (lda1>m || lda2>m ) {
162: for (j=0; j<n; j++) {
163: PetscMemcpy(l->v+j*lda2,mat->v+j*lda1,m*sizeof(PetscScalar));
164: }
165: } else {
166: PetscMemcpy(l->v,mat->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));
167: }
168: (*fact)->factor = 0;
169: MatLUFactor(*fact,0,0,&info);
170: return(0);
171: }
175: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,MatFactorInfo *info,Mat *fact)
176: {
180: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,fact);
181: return(0);
182: }
186: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,MatFactorInfo *factinfo)
187: {
188: #if defined(PETSC_MISSING_LAPACK_POTRF)
190: SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
191: #else
192: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
194: PetscBLASInt n = (PetscBLASInt)A->cmap.n,info;
195:
197: PetscFree(mat->pivots);
198: mat->pivots = 0;
200: if (!A->rmap.n || !A->cmap.n) return(0);
201: LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info);
202: if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
203: A->factor = FACTOR_CHOLESKY;
204: PetscLogFlops((A->cmap.n*A->cmap.n*A->cmap.n)/3);
205: #endif
206: return(0);
207: }
211: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact)
212: {
214: MatFactorInfo info;
217: info.fill = 1.0;
218: MatCholeskyFactor_SeqDense(*fact,0,&info);
219: return(0);
220: }
224: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
225: {
226: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
228: PetscBLASInt m = (PetscBLASInt)A->rmap.n, one = 1,info;
229: PetscScalar *x,*y;
230:
232: VecGetArray(xx,&x);
233: VecGetArray(yy,&y);
234: PetscMemcpy(y,x,A->rmap.n*sizeof(PetscScalar));
235: if (A->factor == FACTOR_LU) {
236: #if defined(PETSC_MISSING_LAPACK_GETRS)
237: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
238: #else
239: LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
240: if (info) SETERRQ(PETSC_ERR_LIB,"GETRS - Bad solve");
241: #endif
242: } else if (A->factor == FACTOR_CHOLESKY){
243: #if defined(PETSC_MISSING_LAPACK_POTRS)
244: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
245: #else
246: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
247: if (info) SETERRQ(PETSC_ERR_LIB,"POTRS Bad solve");
248: #endif
249: }
250: else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
251: VecRestoreArray(xx,&x);
252: VecRestoreArray(yy,&y);
253: PetscLogFlops(2*A->cmap.n*A->cmap.n - A->cmap.n);
254: return(0);
255: }
259: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
260: {
261: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
263: PetscBLASInt m = (PetscBLASInt) A->rmap.n,one = 1,info;
264: PetscScalar *x,*y;
265:
267: VecGetArray(xx,&x);
268: VecGetArray(yy,&y);
269: PetscMemcpy(y,x,A->rmap.n*sizeof(PetscScalar));
270: /* assume if pivots exist then use LU; else Cholesky */
271: if (mat->pivots) {
272: #if defined(PETSC_MISSING_LAPACK_GETRS)
273: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
274: #else
275: LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
276: if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve");
277: #endif
278: } else {
279: #if defined(PETSC_MISSING_LAPACK_POTRS)
280: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
281: #else
282: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
283: if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve");
284: #endif
285: }
286: VecRestoreArray(xx,&x);
287: VecRestoreArray(yy,&y);
288: PetscLogFlops(2*A->cmap.n*A->cmap.n - A->cmap.n);
289: return(0);
290: }
294: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
295: {
296: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
298: PetscBLASInt m = (PetscBLASInt)A->rmap.n,one = 1,info;
299: PetscScalar *x,*y,sone = 1.0;
300: Vec tmp = 0;
301:
303: VecGetArray(xx,&x);
304: VecGetArray(yy,&y);
305: if (!A->rmap.n || !A->cmap.n) return(0);
306: if (yy == zz) {
307: VecDuplicate(yy,&tmp);
308: PetscLogObjectParent(A,tmp);
309: VecCopy(yy,tmp);
310: }
311: PetscMemcpy(y,x,A->rmap.n*sizeof(PetscScalar));
312: /* assume if pivots exist then use LU; else Cholesky */
313: if (mat->pivots) {
314: #if defined(PETSC_MISSING_LAPACK_GETRS)
315: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
316: #else
317: LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
318: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
319: #endif
320: } else {
321: #if defined(PETSC_MISSING_LAPACK_POTRS)
322: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
323: #else
324: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
325: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
326: #endif
327: }
328: if (tmp) {VecAXPY(yy,sone,tmp); VecDestroy(tmp);}
329: else {VecAXPY(yy,sone,zz);}
330: VecRestoreArray(xx,&x);
331: VecRestoreArray(yy,&y);
332: PetscLogFlops(2*A->cmap.n*A->cmap.n);
333: return(0);
334: }
338: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
339: {
340: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
342: PetscBLASInt m = (PetscBLASInt)A->rmap.n,one = 1,info;
343: PetscScalar *x,*y,sone = 1.0;
344: Vec tmp;
345:
347: if (!A->rmap.n || !A->cmap.n) return(0);
348: VecGetArray(xx,&x);
349: VecGetArray(yy,&y);
350: if (yy == zz) {
351: VecDuplicate(yy,&tmp);
352: PetscLogObjectParent(A,tmp);
353: VecCopy(yy,tmp);
354: }
355: PetscMemcpy(y,x,A->rmap.n*sizeof(PetscScalar));
356: /* assume if pivots exist then use LU; else Cholesky */
357: if (mat->pivots) {
358: #if defined(PETSC_MISSING_LAPACK_GETRS)
359: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
360: #else
361: LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
362: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
363: #endif
364: } else {
365: #if defined(PETSC_MISSING_LAPACK_POTRS)
366: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
367: #else
368: LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
369: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
370: #endif
371: }
372: if (tmp) {
373: VecAXPY(yy,sone,tmp);
374: VecDestroy(tmp);
375: } else {
376: VecAXPY(yy,sone,zz);
377: }
378: VecRestoreArray(xx,&x);
379: VecRestoreArray(yy,&y);
380: PetscLogFlops(2*A->cmap.n*A->cmap.n);
381: return(0);
382: }
383: /* ------------------------------------------------------------------*/
386: PetscErrorCode MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
387: {
388: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
389: PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt;
391: PetscInt m = A->rmap.n,i;
392: #if !defined(PETSC_USE_COMPLEX)
393: PetscBLASInt bm = (PetscBLASInt)m, o = 1;
394: #endif
397: if (flag & SOR_ZERO_INITIAL_GUESS) {
398: /* this is a hack fix, should have another version without the second BLASdot */
399: VecSet(xx,zero);
400: }
401: VecGetArray(xx,&x);
402: VecGetArray(bb,&b);
403: its = its*lits;
404: if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
405: while (its--) {
406: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
407: for (i=0; i<m; i++) {
408: #if defined(PETSC_USE_COMPLEX)
409: /* cannot use BLAS dot for complex because compiler/linker is
410: not happy about returning a double complex */
411: PetscInt _i;
412: PetscScalar sum = b[i];
413: for (_i=0; _i<m; _i++) {
414: sum -= PetscConj(v[i+_i*m])*x[_i];
415: }
416: xt = sum;
417: #else
418: xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
419: #endif
420: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
421: }
422: }
423: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
424: for (i=m-1; i>=0; i--) {
425: #if defined(PETSC_USE_COMPLEX)
426: /* cannot use BLAS dot for complex because compiler/linker is
427: not happy about returning a double complex */
428: PetscInt _i;
429: PetscScalar sum = b[i];
430: for (_i=0; _i<m; _i++) {
431: sum -= PetscConj(v[i+_i*m])*x[_i];
432: }
433: xt = sum;
434: #else
435: xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
436: #endif
437: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
438: }
439: }
440: }
441: VecRestoreArray(bb,&b);
442: VecRestoreArray(xx,&x);
443: return(0);
444: }
446: /* -----------------------------------------------------------------*/
449: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
450: {
451: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
452: PetscScalar *v = mat->v,*x,*y;
454: PetscBLASInt m = (PetscBLASInt)A->rmap.n, n = (PetscBLASInt)A->cmap.n,_One=1;
455: PetscScalar _DOne=1.0,_DZero=0.0;
458: if (!A->rmap.n || !A->cmap.n) return(0);
459: VecGetArray(xx,&x);
460: VecGetArray(yy,&y);
461: BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One);
462: VecRestoreArray(xx,&x);
463: VecRestoreArray(yy,&y);
464: PetscLogFlops(2*A->rmap.n*A->cmap.n - A->cmap.n);
465: return(0);
466: }
470: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
471: {
472: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
473: PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
475: PetscBLASInt m = (PetscBLASInt)A->rmap.n, n = (PetscBLASInt)A->cmap.n, _One=1;
478: if (!A->rmap.n || !A->cmap.n) return(0);
479: VecGetArray(xx,&x);
480: VecGetArray(yy,&y);
481: BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One);
482: VecRestoreArray(xx,&x);
483: VecRestoreArray(yy,&y);
484: PetscLogFlops(2*A->rmap.n*A->cmap.n - A->rmap.n);
485: return(0);
486: }
490: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
491: {
492: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
493: PetscScalar *v = mat->v,*x,*y,_DOne=1.0;
495: PetscBLASInt m = (PetscBLASInt)A->rmap.n, n = (PetscBLASInt)A->cmap.n, _One=1;
498: if (!A->rmap.n || !A->cmap.n) return(0);
499: if (zz != yy) {VecCopy(zz,yy);}
500: VecGetArray(xx,&x);
501: VecGetArray(yy,&y);
502: BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
503: VecRestoreArray(xx,&x);
504: VecRestoreArray(yy,&y);
505: PetscLogFlops(2*A->rmap.n*A->cmap.n);
506: return(0);
507: }
511: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
512: {
513: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
514: PetscScalar *v = mat->v,*x,*y;
516: PetscBLASInt m = (PetscBLASInt)A->rmap.n, n = (PetscBLASInt)A->cmap.n, _One=1;
517: PetscScalar _DOne=1.0;
520: if (!A->rmap.n || !A->cmap.n) return(0);
521: if (zz != yy) {VecCopy(zz,yy);}
522: VecGetArray(xx,&x);
523: VecGetArray(yy,&y);
524: BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
525: VecRestoreArray(xx,&x);
526: VecRestoreArray(yy,&y);
527: PetscLogFlops(2*A->rmap.n*A->cmap.n);
528: return(0);
529: }
531: /* -----------------------------------------------------------------*/
534: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
535: {
536: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
537: PetscScalar *v;
539: PetscInt i;
540:
542: *ncols = A->cmap.n;
543: if (cols) {
544: PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),cols);
545: for (i=0; i<A->cmap.n; i++) (*cols)[i] = i;
546: }
547: if (vals) {
548: PetscMalloc((A->cmap.n+1)*sizeof(PetscScalar),vals);
549: v = mat->v + row;
550: for (i=0; i<A->cmap.n; i++) {(*vals)[i] = *v; v += mat->lda;}
551: }
552: return(0);
553: }
557: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
558: {
561: if (cols) {PetscFree(*cols);}
562: if (vals) {PetscFree(*vals); }
563: return(0);
564: }
565: /* ----------------------------------------------------------------*/
568: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
569: {
570: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
571: PetscInt i,j,idx=0;
572:
574: if (!mat->roworiented) {
575: if (addv == INSERT_VALUES) {
576: for (j=0; j<n; j++) {
577: if (indexn[j] < 0) {idx += m; continue;}
578: #if defined(PETSC_USE_DEBUG)
579: if (indexn[j] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap.n-1);
580: #endif
581: for (i=0; i<m; i++) {
582: if (indexm[i] < 0) {idx++; continue;}
583: #if defined(PETSC_USE_DEBUG)
584: if (indexm[i] >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap.n-1);
585: #endif
586: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
587: }
588: }
589: } else {
590: for (j=0; j<n; j++) {
591: if (indexn[j] < 0) {idx += m; continue;}
592: #if defined(PETSC_USE_DEBUG)
593: if (indexn[j] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap.n-1);
594: #endif
595: for (i=0; i<m; i++) {
596: if (indexm[i] < 0) {idx++; continue;}
597: #if defined(PETSC_USE_DEBUG)
598: if (indexm[i] >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap.n-1);
599: #endif
600: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
601: }
602: }
603: }
604: } else {
605: if (addv == INSERT_VALUES) {
606: for (i=0; i<m; i++) {
607: if (indexm[i] < 0) { idx += n; continue;}
608: #if defined(PETSC_USE_DEBUG)
609: if (indexm[i] >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap.n-1);
610: #endif
611: for (j=0; j<n; j++) {
612: if (indexn[j] < 0) { idx++; continue;}
613: #if defined(PETSC_USE_DEBUG)
614: if (indexn[j] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap.n-1);
615: #endif
616: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
617: }
618: }
619: } else {
620: for (i=0; i<m; i++) {
621: if (indexm[i] < 0) { idx += n; continue;}
622: #if defined(PETSC_USE_DEBUG)
623: if (indexm[i] >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap.n-1);
624: #endif
625: for (j=0; j<n; j++) {
626: if (indexn[j] < 0) { idx++; continue;}
627: #if defined(PETSC_USE_DEBUG)
628: if (indexn[j] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap.n-1);
629: #endif
630: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
631: }
632: }
633: }
634: }
635: return(0);
636: }
640: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
641: {
642: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
643: PetscInt i,j;
644: PetscScalar *vpt = v;
647: /* row-oriented output */
648: for (i=0; i<m; i++) {
649: for (j=0; j<n; j++) {
650: *vpt++ = mat->v[indexn[j]*mat->lda + indexm[i]];
651: }
652: }
653: return(0);
654: }
656: /* -----------------------------------------------------------------*/
658: #include petscsys.h
662: PetscErrorCode MatLoad_SeqDense(PetscViewer viewer, MatType type,Mat *A)
663: {
664: Mat_SeqDense *a;
665: Mat B;
667: PetscInt *scols,i,j,nz,header[4];
668: int fd;
669: PetscMPIInt size;
670: PetscInt *rowlengths = 0,M,N,*cols;
671: PetscScalar *vals,*svals,*v,*w;
672: MPI_Comm comm = ((PetscObject)viewer)->comm;
675: MPI_Comm_size(comm,&size);
676: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor");
677: PetscViewerBinaryGetDescriptor(viewer,&fd);
678: PetscBinaryRead(fd,header,4,PETSC_INT);
679: if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
680: M = header[1]; N = header[2]; nz = header[3];
682: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
683: MatCreate(comm,A);
684: MatSetSizes(*A,M,N,M,N);
685: MatSetType(*A,type);
686: MatSeqDenseSetPreallocation(*A,PETSC_NULL);
687: B = *A;
688: a = (Mat_SeqDense*)B->data;
689: v = a->v;
690: /* Allocate some temp space to read in the values and then flip them
691: from row major to column major */
692: PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);
693: /* read in nonzero values */
694: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
695: /* now flip the values and store them in the matrix*/
696: for (j=0; j<N; j++) {
697: for (i=0; i<M; i++) {
698: *v++ =w[i*N+j];
699: }
700: }
701: PetscFree(w);
702: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
703: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
704: } else {
705: /* read row lengths */
706: PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);
707: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
709: /* create our matrix */
710: MatCreate(comm,A);
711: MatSetSizes(*A,M,N,M,N);
712: MatSetType(*A,type);
713: MatSeqDenseSetPreallocation(*A,PETSC_NULL);
714: B = *A;
715: a = (Mat_SeqDense*)B->data;
716: v = a->v;
718: /* read column indices and nonzeros */
719: PetscMalloc((nz+1)*sizeof(PetscInt),&scols);
720: cols = scols;
721: PetscBinaryRead(fd,cols,nz,PETSC_INT);
722: PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);
723: vals = svals;
724: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
726: /* insert into matrix */
727: for (i=0; i<M; i++) {
728: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
729: svals += rowlengths[i]; scols += rowlengths[i];
730: }
731: PetscFree(vals);
732: PetscFree(cols);
733: PetscFree(rowlengths);
735: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
736: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
737: }
738: return(0);
739: }
741: #include petscsys.h
745: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
746: {
747: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
748: PetscErrorCode ierr;
749: PetscInt i,j;
750: const char *name;
751: PetscScalar *v;
752: PetscViewerFormat format;
753: #if defined(PETSC_USE_COMPLEX)
754: PetscTruth allreal = PETSC_TRUE;
755: #endif
758: PetscObjectGetName((PetscObject)A,&name);
759: PetscViewerGetFormat(viewer,&format);
760: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
761: return(0); /* do nothing for now */
762: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
763: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
764: for (i=0; i<A->rmap.n; i++) {
765: v = a->v + i;
766: PetscViewerASCIIPrintf(viewer,"row %D:",i);
767: for (j=0; j<A->cmap.n; j++) {
768: #if defined(PETSC_USE_COMPLEX)
769: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
770: PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));
771: } else if (PetscRealPart(*v)) {
772: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));
773: }
774: #else
775: if (*v) {
776: PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);
777: }
778: #endif
779: v += a->lda;
780: }
781: PetscViewerASCIIPrintf(viewer,"\n");
782: }
783: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
784: } else {
785: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
786: #if defined(PETSC_USE_COMPLEX)
787: /* determine if matrix has all real values */
788: v = a->v;
789: for (i=0; i<A->rmap.n*A->cmap.n; i++) {
790: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;}
791: }
792: #endif
793: if (format == PETSC_VIEWER_ASCII_MATLAB) {
794: PetscObjectGetName((PetscObject)A,&name);
795: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap.n,A->cmap.n);
796: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap.n,A->cmap.n);
797: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
798: }
800: for (i=0; i<A->rmap.n; i++) {
801: v = a->v + i;
802: for (j=0; j<A->cmap.n; j++) {
803: #if defined(PETSC_USE_COMPLEX)
804: if (allreal) {
805: PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));
806: } else {
807: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));
808: }
809: #else
810: PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);
811: #endif
812: v += a->lda;
813: }
814: PetscViewerASCIIPrintf(viewer,"\n");
815: }
816: if (format == PETSC_VIEWER_ASCII_MATLAB) {
817: PetscViewerASCIIPrintf(viewer,"];\n");
818: }
819: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
820: }
821: PetscViewerFlush(viewer);
822: return(0);
823: }
827: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
828: {
829: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
830: PetscErrorCode ierr;
831: int fd;
832: PetscInt ict,j,n = A->cmap.n,m = A->rmap.n,i,*col_lens,nz = m*n;
833: PetscScalar *v,*anonz,*vals;
834: PetscViewerFormat format;
835:
837: PetscViewerBinaryGetDescriptor(viewer,&fd);
839: PetscViewerGetFormat(viewer,&format);
840: if (format == PETSC_VIEWER_BINARY_NATIVE) {
841: /* store the matrix as a dense matrix */
842: PetscMalloc(4*sizeof(PetscInt),&col_lens);
843: col_lens[0] = MAT_FILE_COOKIE;
844: col_lens[1] = m;
845: col_lens[2] = n;
846: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
847: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
848: PetscFree(col_lens);
850: /* write out matrix, by rows */
851: PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);
852: v = a->v;
853: for (i=0; i<m; i++) {
854: for (j=0; j<n; j++) {
855: vals[i + j*m] = *v++;
856: }
857: }
858: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
859: PetscFree(vals);
860: } else {
861: PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);
862: col_lens[0] = MAT_FILE_COOKIE;
863: col_lens[1] = m;
864: col_lens[2] = n;
865: col_lens[3] = nz;
867: /* store lengths of each row and write (including header) to file */
868: for (i=0; i<m; i++) col_lens[4+i] = n;
869: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
871: /* Possibly should write in smaller increments, not whole matrix at once? */
872: /* store column indices (zero start index) */
873: ict = 0;
874: for (i=0; i<m; i++) {
875: for (j=0; j<n; j++) col_lens[ict++] = j;
876: }
877: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
878: PetscFree(col_lens);
880: /* store nonzero values */
881: PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);
882: ict = 0;
883: for (i=0; i<m; i++) {
884: v = a->v + i;
885: for (j=0; j<n; j++) {
886: anonz[ict++] = *v; v += a->lda;
887: }
888: }
889: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
890: PetscFree(anonz);
891: }
892: return(0);
893: }
897: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
898: {
899: Mat A = (Mat) Aa;
900: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
901: PetscErrorCode ierr;
902: PetscInt m = A->rmap.n,n = A->cmap.n,color,i,j;
903: PetscScalar *v = a->v;
904: PetscViewer viewer;
905: PetscDraw popup;
906: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
907: PetscViewerFormat format;
911: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
912: PetscViewerGetFormat(viewer,&format);
913: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
915: /* Loop over matrix elements drawing boxes */
916: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
917: /* Blue for negative and Red for positive */
918: color = PETSC_DRAW_BLUE;
919: for(j = 0; j < n; j++) {
920: x_l = j;
921: x_r = x_l + 1.0;
922: for(i = 0; i < m; i++) {
923: y_l = m - i - 1.0;
924: y_r = y_l + 1.0;
925: #if defined(PETSC_USE_COMPLEX)
926: if (PetscRealPart(v[j*m+i]) > 0.) {
927: color = PETSC_DRAW_RED;
928: } else if (PetscRealPart(v[j*m+i]) < 0.) {
929: color = PETSC_DRAW_BLUE;
930: } else {
931: continue;
932: }
933: #else
934: if (v[j*m+i] > 0.) {
935: color = PETSC_DRAW_RED;
936: } else if (v[j*m+i] < 0.) {
937: color = PETSC_DRAW_BLUE;
938: } else {
939: continue;
940: }
941: #endif
942: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
943: }
944: }
945: } else {
946: /* use contour shading to indicate magnitude of values */
947: /* first determine max of all nonzero values */
948: for(i = 0; i < m*n; i++) {
949: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
950: }
951: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
952: PetscDrawGetPopup(draw,&popup);
953: if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
954: for(j = 0; j < n; j++) {
955: x_l = j;
956: x_r = x_l + 1.0;
957: for(i = 0; i < m; i++) {
958: y_l = m - i - 1.0;
959: y_r = y_l + 1.0;
960: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
961: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
962: }
963: }
964: }
965: return(0);
966: }
970: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
971: {
972: PetscDraw draw;
973: PetscTruth isnull;
974: PetscReal xr,yr,xl,yl,h,w;
978: PetscViewerDrawGetDraw(viewer,0,&draw);
979: PetscDrawIsNull(draw,&isnull);
980: if (isnull) return(0);
982: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
983: xr = A->cmap.n; yr = A->rmap.n; h = yr/10.0; w = xr/10.0;
984: xr += w; yr += h; xl = -w; yl = -h;
985: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
986: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
987: PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
988: return(0);
989: }
993: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
994: {
996: PetscTruth iascii,isbinary,isdraw;
999: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1000: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
1001: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
1003: if (iascii) {
1004: MatView_SeqDense_ASCII(A,viewer);
1005: } else if (isbinary) {
1006: MatView_SeqDense_Binary(A,viewer);
1007: } else if (isdraw) {
1008: MatView_SeqDense_Draw(A,viewer);
1009: } else {
1010: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name);
1011: }
1012: return(0);
1013: }
1017: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1018: {
1019: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1023: #if defined(PETSC_USE_LOG)
1024: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap.n,mat->cmap.n);
1025: #endif
1026: PetscFree(l->pivots);
1027: if (!l->user_alloc) {PetscFree(l->v);}
1028: PetscFree(l);
1030: PetscObjectChangeTypeName((PetscObject)mat,0);
1031: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);
1032: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);
1033: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);
1034: PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);
1035: return(0);
1036: }
1040: PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout)
1041: {
1042: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1044: PetscInt k,j,m,n,M;
1045: PetscScalar *v,tmp;
1048: v = mat->v; m = A->rmap.n; M = mat->lda; n = A->cmap.n;
1049: if (!matout) { /* in place transpose */
1050: if (m != n) {
1051: SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1052: } else {
1053: for (j=0; j<m; j++) {
1054: for (k=0; k<j; k++) {
1055: tmp = v[j + k*M];
1056: v[j + k*M] = v[k + j*M];
1057: v[k + j*M] = tmp;
1058: }
1059: }
1060: }
1061: } else { /* out-of-place transpose */
1062: Mat tmat;
1063: Mat_SeqDense *tmatd;
1064: PetscScalar *v2;
1066: MatCreate(A->comm,&tmat);
1067: MatSetSizes(tmat,A->cmap.n,A->rmap.n,A->cmap.n,A->rmap.n);
1068: MatSetType(tmat,A->type_name);
1069: MatSeqDenseSetPreallocation(tmat,PETSC_NULL);
1070: tmatd = (Mat_SeqDense*)tmat->data;
1071: v = mat->v; v2 = tmatd->v;
1072: for (j=0; j<n; j++) {
1073: for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
1074: }
1075: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1076: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1077: *matout = tmat;
1078: }
1079: return(0);
1080: }
1084: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg)
1085: {
1086: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1087: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1088: PetscInt i,j;
1089: PetscScalar *v1 = mat1->v,*v2 = mat2->v;
1092: if (A1->rmap.n != A2->rmap.n) {*flg = PETSC_FALSE; return(0);}
1093: if (A1->cmap.n != A2->cmap.n) {*flg = PETSC_FALSE; return(0);}
1094: for (i=0; i<A1->rmap.n; i++) {
1095: v1 = mat1->v+i; v2 = mat2->v+i;
1096: for (j=0; j<A1->cmap.n; j++) {
1097: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1098: v1 += mat1->lda; v2 += mat2->lda;
1099: }
1100: }
1101: *flg = PETSC_TRUE;
1102: return(0);
1103: }
1107: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1108: {
1109: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1111: PetscInt i,n,len;
1112: PetscScalar *x,zero = 0.0;
1115: VecSet(v,zero);
1116: VecGetSize(v,&n);
1117: VecGetArray(v,&x);
1118: len = PetscMin(A->rmap.n,A->cmap.n);
1119: if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1120: for (i=0; i<len; i++) {
1121: x[i] = mat->v[i*mat->lda + i];
1122: }
1123: VecRestoreArray(v,&x);
1124: return(0);
1125: }
1129: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1130: {
1131: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1132: PetscScalar *l,*r,x,*v;
1134: PetscInt i,j,m = A->rmap.n,n = A->cmap.n;
1137: if (ll) {
1138: VecGetSize(ll,&m);
1139: VecGetArray(ll,&l);
1140: if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1141: for (i=0; i<m; i++) {
1142: x = l[i];
1143: v = mat->v + i;
1144: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1145: }
1146: VecRestoreArray(ll,&l);
1147: PetscLogFlops(n*m);
1148: }
1149: if (rr) {
1150: VecGetSize(rr,&n);
1151: VecGetArray(rr,&r);
1152: if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1153: for (i=0; i<n; i++) {
1154: x = r[i];
1155: v = mat->v + i*m;
1156: for (j=0; j<m; j++) { (*v++) *= x;}
1157: }
1158: VecRestoreArray(rr,&r);
1159: PetscLogFlops(n*m);
1160: }
1161: return(0);
1162: }
1166: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1167: {
1168: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1169: PetscScalar *v = mat->v;
1170: PetscReal sum = 0.0;
1171: PetscInt lda=mat->lda,m=A->rmap.n,i,j;
1175: if (type == NORM_FROBENIUS) {
1176: if (lda>m) {
1177: for (j=0; j<A->cmap.n; j++) {
1178: v = mat->v+j*lda;
1179: for (i=0; i<m; i++) {
1180: #if defined(PETSC_USE_COMPLEX)
1181: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1182: #else
1183: sum += (*v)*(*v); v++;
1184: #endif
1185: }
1186: }
1187: } else {
1188: for (i=0; i<A->cmap.n*A->rmap.n; i++) {
1189: #if defined(PETSC_USE_COMPLEX)
1190: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1191: #else
1192: sum += (*v)*(*v); v++;
1193: #endif
1194: }
1195: }
1196: *nrm = sqrt(sum);
1197: PetscLogFlops(2*A->cmap.n*A->rmap.n);
1198: } else if (type == NORM_1) {
1199: *nrm = 0.0;
1200: for (j=0; j<A->cmap.n; j++) {
1201: v = mat->v + j*mat->lda;
1202: sum = 0.0;
1203: for (i=0; i<A->rmap.n; i++) {
1204: sum += PetscAbsScalar(*v); v++;
1205: }
1206: if (sum > *nrm) *nrm = sum;
1207: }
1208: PetscLogFlops(A->cmap.n*A->rmap.n);
1209: } else if (type == NORM_INFINITY) {
1210: *nrm = 0.0;
1211: for (j=0; j<A->rmap.n; j++) {
1212: v = mat->v + j;
1213: sum = 0.0;
1214: for (i=0; i<A->cmap.n; i++) {
1215: sum += PetscAbsScalar(*v); v += mat->lda;
1216: }
1217: if (sum > *nrm) *nrm = sum;
1218: }
1219: PetscLogFlops(A->cmap.n*A->rmap.n);
1220: } else {
1221: SETERRQ(PETSC_ERR_SUP,"No two norm");
1222: }
1223: return(0);
1224: }
1228: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op)
1229: {
1230: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1232:
1234: switch (op) {
1235: case MAT_ROW_ORIENTED:
1236: aij->roworiented = PETSC_TRUE;
1237: break;
1238: case MAT_COLUMN_ORIENTED:
1239: aij->roworiented = PETSC_FALSE;
1240: break;
1241: case MAT_ROWS_SORTED:
1242: case MAT_ROWS_UNSORTED:
1243: case MAT_COLUMNS_SORTED:
1244: case MAT_COLUMNS_UNSORTED:
1245: case MAT_NO_NEW_NONZERO_LOCATIONS:
1246: case MAT_YES_NEW_NONZERO_LOCATIONS:
1247: case MAT_NEW_NONZERO_LOCATION_ERR:
1248: case MAT_NO_NEW_DIAGONALS:
1249: case MAT_YES_NEW_DIAGONALS:
1250: case MAT_IGNORE_OFF_PROC_ENTRIES:
1251: case MAT_USE_HASH_TABLE:
1252: case MAT_SYMMETRIC:
1253: case MAT_STRUCTURALLY_SYMMETRIC:
1254: case MAT_NOT_SYMMETRIC:
1255: case MAT_NOT_STRUCTURALLY_SYMMETRIC:
1256: case MAT_HERMITIAN:
1257: case MAT_NOT_HERMITIAN:
1258: case MAT_SYMMETRY_ETERNAL:
1259: case MAT_NOT_SYMMETRY_ETERNAL:
1260: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1261: break;
1262: default:
1263: SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op);
1264: }
1265: return(0);
1266: }
1270: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1271: {
1272: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1274: PetscInt lda=l->lda,m=A->rmap.n,j;
1277: if (lda>m) {
1278: for (j=0; j<A->cmap.n; j++) {
1279: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1280: }
1281: } else {
1282: PetscMemzero(l->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));
1283: }
1284: return(0);
1285: }
1289: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscI