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 PetscInt rows[],PetscScalar diag)
1290: {
1291: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1292: PetscInt n = A->cmap.n,i,j;
1293: PetscScalar *slot;
1296: for (i=0; i<N; i++) {
1297: slot = l->v + rows[i];
1298: for (j=0; j<n; j++) { *slot = 0.0; slot += n;}
1299: }
1300: if (diag != 0.0) {
1301: for (i=0; i<N; i++) {
1302: slot = l->v + (n+1)*rows[i];
1303: *slot = diag;
1304: }
1305: }
1306: return(0);
1307: }
1311: PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[])
1312: {
1313: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1316: if (mat->lda != A->rmap.n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows");
1317: *array = mat->v;
1318: return(0);
1319: }
1323: PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1324: {
1326: *array = 0; /* user cannot accidently use the array later */
1327: return(0);
1328: }
1332: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1333: {
1334: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1336: PetscInt i,j,*irow,*icol,nrows,ncols;
1337: PetscScalar *av,*bv,*v = mat->v;
1338: Mat newmat;
1341: ISGetIndices(isrow,&irow);
1342: ISGetIndices(iscol,&icol);
1343: ISGetLocalSize(isrow,&nrows);
1344: ISGetLocalSize(iscol,&ncols);
1345:
1346: /* Check submatrixcall */
1347: if (scall == MAT_REUSE_MATRIX) {
1348: PetscInt n_cols,n_rows;
1349: MatGetSize(*B,&n_rows,&n_cols);
1350: if (n_rows != nrows || n_cols != ncols) {
1351: /* resize the result result matrix to match number of requested rows/columns */
1352: MatSetSizes(*B,nrows,nrows,nrows,nrows);
1353: }
1354: newmat = *B;
1355: } else {
1356: /* Create and fill new matrix */
1357: MatCreate(A->comm,&newmat);
1358: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1359: MatSetType(newmat,A->type_name);
1360: MatSeqDenseSetPreallocation(newmat,PETSC_NULL);
1361: }
1363: /* Now extract the data pointers and do the copy,column at a time */
1364: bv = ((Mat_SeqDense*)newmat->data)->v;
1365:
1366: for (i=0; i<ncols; i++) {
1367: av = v + mat->lda*icol[i];
1368: for (j=0; j<nrows; j++) {
1369: *bv++ = av[irow[j]];
1370: }
1371: }
1373: /* Assemble the matrices so that the correct flags are set */
1374: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1375: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1377: /* Free work space */
1378: ISRestoreIndices(isrow,&irow);
1379: ISRestoreIndices(iscol,&icol);
1380: *B = newmat;
1381: return(0);
1382: }
1386: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1387: {
1389: PetscInt i;
1392: if (scall == MAT_INITIAL_MATRIX) {
1393: PetscMalloc((n+1)*sizeof(Mat),B);
1394: }
1396: for (i=0; i<n; i++) {
1397: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1398: }
1399: return(0);
1400: }
1404: PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1405: {
1407: return(0);
1408: }
1412: PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1413: {
1415: return(0);
1416: }
1420: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1421: {
1422: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data;
1424: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap.n,n=A->cmap.n, j;
1427: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1428: if (A->ops->copy != B->ops->copy) {
1429: MatCopy_Basic(A,B,str);
1430: return(0);
1431: }
1432: if (m != B->rmap.n || n != B->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1433: if (lda1>m || lda2>m) {
1434: for (j=0; j<n; j++) {
1435: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1436: }
1437: } else {
1438: PetscMemcpy(b->v,a->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));
1439: }
1440: return(0);
1441: }
1445: PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A)
1446: {
1450: MatSeqDenseSetPreallocation(A,0);
1451: return(0);
1452: }
1456: PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N)
1457: {
1458: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1461: /* this will not be called before lda, Mmax, and Nmax have been set */
1462: m = PetscMax(m,M);
1463: n = PetscMax(n,N);
1464: if (m > a->Mmax) SETERRQ2(PETSC_ERR_SUP,"Cannot yet resize number rows of dense matrix larger then its initial size %d, requested %d",a->lda,(int)m);
1465: if (n > a->Nmax) SETERRQ2(PETSC_ERR_SUP,"Cannot yet resize number columns of dense matrix larger then its initial size %d, requested %d",a->Nmax,(int)n);
1466: A->rmap.n = A->rmap.n = m;
1467: A->cmap.n = A->cmap.N = n;
1468: if (a->changelda) a->lda = m;
1469: PetscMemzero(a->v,m*n*sizeof(PetscScalar));
1470: return(0);
1471: }
1473: /* ----------------------------------------------------------------*/
1476: PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1477: {
1481: if (scall == MAT_INITIAL_MATRIX){
1482: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1483: }
1484: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1485: return(0);
1486: }
1490: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1491: {
1493: PetscInt m=A->rmap.n,n=B->cmap.n;
1494: Mat Cmat;
1497: if (A->cmap.n != B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->cmap.n %d != B->rmap.n %d\n",A->cmap.n,B->rmap.n);
1498: MatCreate(PETSC_COMM_SELF,&Cmat);
1499: MatSetSizes(Cmat,m,n,m,n);
1500: MatSetType(Cmat,MATSEQDENSE);
1501: MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);
1502: Cmat->assembled = PETSC_TRUE;
1503: *C = Cmat;
1504: return(0);
1505: }
1509: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1510: {
1511: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1512: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1513: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1514: PetscBLASInt m=(PetscBLASInt)A->rmap.n,n=(PetscBLASInt)B->cmap.n,k=(PetscBLASInt)A->cmap.n;
1515: PetscScalar _DOne=1.0,_DZero=0.0;
1518: BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda);
1519: return(0);
1520: }
1524: PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1525: {
1529: if (scall == MAT_INITIAL_MATRIX){
1530: MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);
1531: }
1532: MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);
1533: return(0);
1534: }
1538: PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1539: {
1541: PetscInt m=A->cmap.n,n=B->cmap.n;
1542: Mat Cmat;
1545: if (A->rmap.n != B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->rmap.n %d != B->rmap.n %d\n",A->rmap.n,B->rmap.n);
1546: MatCreate(PETSC_COMM_SELF,&Cmat);
1547: MatSetSizes(Cmat,m,n,m,n);
1548: MatSetType(Cmat,MATSEQDENSE);
1549: MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);
1550: Cmat->assembled = PETSC_TRUE;
1551: *C = Cmat;
1552: return(0);
1553: }
1557: PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1558: {
1559: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1560: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1561: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1562: PetscBLASInt m=(PetscBLASInt)A->cmap.n,n=(PetscBLASInt)B->cmap.n,k=(PetscBLASInt)A->rmap.n;
1563: PetscScalar _DOne=1.0,_DZero=0.0;
1566: BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda);
1567: return(0);
1568: }
1569: /* -------------------------------------------------------------------*/
1570: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
1571: MatGetRow_SeqDense,
1572: MatRestoreRow_SeqDense,
1573: MatMult_SeqDense,
1574: /* 4*/ MatMultAdd_SeqDense,
1575: MatMultTranspose_SeqDense,
1576: MatMultTransposeAdd_SeqDense,
1577: MatSolve_SeqDense,
1578: MatSolveAdd_SeqDense,
1579: MatSolveTranspose_SeqDense,
1580: /*10*/ MatSolveTransposeAdd_SeqDense,
1581: MatLUFactor_SeqDense,
1582: MatCholeskyFactor_SeqDense,
1583: MatRelax_SeqDense,
1584: MatTranspose_SeqDense,
1585: /*15*/ MatGetInfo_SeqDense,
1586: MatEqual_SeqDense,
1587: MatGetDiagonal_SeqDense,
1588: MatDiagonalScale_SeqDense,
1589: MatNorm_SeqDense,
1590: /*20*/ MatAssemblyBegin_SeqDense,
1591: MatAssemblyEnd_SeqDense,
1592: 0,
1593: MatSetOption_SeqDense,
1594: MatZeroEntries_SeqDense,
1595: /*25*/ MatZeroRows_SeqDense,
1596: MatLUFactorSymbolic_SeqDense,
1597: MatLUFactorNumeric_SeqDense,
1598: MatCholeskyFactorSymbolic_SeqDense,
1599: MatCholeskyFactorNumeric_SeqDense,
1600: /*30*/ MatSetUpPreallocation_SeqDense,
1601: 0,
1602: 0,
1603: MatGetArray_SeqDense,
1604: MatRestoreArray_SeqDense,
1605: /*35*/ MatDuplicate_SeqDense,
1606: 0,
1607: 0,
1608: 0,
1609: 0,
1610: /*40*/ MatAXPY_SeqDense,
1611: MatGetSubMatrices_SeqDense,
1612: 0,
1613: MatGetValues_SeqDense,
1614: MatCopy_SeqDense,
1615: /*45*/ 0,
1616: MatScale_SeqDense,
1617: 0,
1618: 0,
1619: 0,
1620: /*50*/ 0,
1621: 0,
1622: 0,
1623: 0,
1624: 0,
1625: /*55*/ 0,
1626: 0,
1627: 0,
1628: 0,
1629: 0,
1630: /*60*/ 0,
1631: MatDestroy_SeqDense,
1632: MatView_SeqDense,
1633: 0,
1634: 0,
1635: /*65*/ 0,
1636: 0,
1637: 0,
1638: 0,
1639: 0,
1640: /*70*/ 0,
1641: 0,
1642: 0,
1643: 0,
1644: 0,
1645: /*75*/ 0,
1646: 0,
1647: 0,
1648: 0,
1649: 0,
1650: /*80*/ 0,
1651: 0,
1652: 0,
1653: 0,
1654: /*84*/ MatLoad_SeqDense,
1655: 0,
1656: 0,
1657: 0,
1658: 0,
1659: 0,
1660: /*90*/ MatMatMult_SeqDense_SeqDense,
1661: MatMatMultSymbolic_SeqDense_SeqDense,
1662: MatMatMultNumeric_SeqDense_SeqDense,
1663: 0,
1664: 0,
1665: /*95*/ 0,
1666: MatMatMultTranspose_SeqDense_SeqDense,
1667: MatMatMultTransposeSymbolic_SeqDense_SeqDense,
1668: MatMatMultTransposeNumeric_SeqDense_SeqDense,
1669: 0,
1670: /*100*/0,
1671: 0,
1672: 0,
1673: 0,
1674: MatSetSizes_SeqDense};
1678: /*@C
1679: MatCreateSeqDense - Creates a sequential dense matrix that
1680: is stored in column major order (the usual Fortran 77 manner). Many
1681: of the matrix operations use the BLAS and LAPACK routines.
1683: Collective on MPI_Comm
1685: Input Parameters:
1686: + comm - MPI communicator, set to PETSC_COMM_SELF
1687: . m - number of rows
1688: . n - number of columns
1689: - data - optional location of matrix data. Set data=PETSC_NULL for PETSc
1690: to control all matrix memory allocation.
1692: Output Parameter:
1693: . A - the matrix
1695: Notes:
1696: The data input variable is intended primarily for Fortran programmers
1697: who wish to allocate their own matrix memory space. Most users should
1698: set data=PETSC_NULL.
1700: Level: intermediate
1702: .keywords: dense, matrix, LAPACK, BLAS
1704: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
1705: @*/
1706: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
1707: {
1711: MatCreate(comm,A);
1712: MatSetSizes(*A,m,n,m,n);
1713: MatSetType(*A,MATSEQDENSE);
1714: MatSeqDenseSetPreallocation(*A,data);
1715: return(0);
1716: }
1720: /*@C
1721: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
1723: Collective on MPI_Comm
1725: Input Parameters:
1726: + A - the matrix
1727: - data - the array (or PETSC_NULL)
1729: Notes:
1730: The data input variable is intended primarily for Fortran programmers
1731: who wish to allocate their own matrix memory space. Most users should
1732: need not call this routine.
1734: Level: intermediate
1736: .keywords: dense, matrix, LAPACK, BLAS
1738: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
1739: @*/
1740: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
1741: {
1742: PetscErrorCode ierr,(*f)(Mat,PetscScalar[]);
1745: PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);
1746: if (f) {
1747: (*f)(B,data);
1748: }
1749: return(0);
1750: }
1755: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
1756: {
1757: Mat_SeqDense *b;
1761: B->preallocated = PETSC_TRUE;
1762: b = (Mat_SeqDense*)B->data;
1763: if (!data) {
1764: PetscMalloc((b->lda*b->Nmax+1)*sizeof(PetscScalar),&b->v);
1765: PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));
1766: b->user_alloc = PETSC_FALSE;
1767: PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));
1768: } else { /* user-allocated storage */
1769: b->v = data;
1770: b->user_alloc = PETSC_TRUE;
1771: }
1772: return(0);
1773: }
1778: /*@C
1779: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
1781: Input parameter:
1782: + A - the matrix
1783: - lda - the leading dimension
1785: Notes:
1786: This routine is to be used in conjunction with MatSeqDenseSetPreallocation;
1787: it asserts that the preallocation has a leading dimension (the LDA parameter
1788: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
1790: Level: intermediate
1792: .keywords: dense, matrix, LAPACK, BLAS
1794: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
1795: @*/
1796: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
1797: {
1798: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1801: if (lda < B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap.n);
1802: b->lda = lda;
1803: b->changelda = PETSC_FALSE;
1804: b->Mmax = PetscMax(b->Mmax,lda);
1805: return(0);
1806: }
1808: /*MC
1809: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
1811: Options Database Keys:
1812: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
1814: Level: beginner
1816: .seealso: MatCreateSeqDense()
1818: M*/
1823: PetscErrorCode MatCreate_SeqDense(Mat B)
1824: {
1825: Mat_SeqDense *b;
1827: PetscMPIInt size;
1830: MPI_Comm_size(B->comm,&size);
1831: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1833: B->rmap.bs = B->cmap.bs = 1;
1834: PetscMapInitialize(B->comm,&B->rmap);
1835: PetscMapInitialize(B->comm,&B->cmap);
1837: PetscNew(Mat_SeqDense,&b);
1838: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1839: B->factor = 0;
1840: B->mapping = 0;
1841: PetscLogObjectMemory(B,sizeof(struct _p_Mat));
1842: B->data = (void*)b;
1845: b->pivots = 0;
1846: b->roworiented = PETSC_TRUE;
1847: b->v = 0;
1848: b->lda = B->rmap.n;
1849: b->changelda = PETSC_FALSE;
1850: b->Mmax = B->rmap.n;
1851: b->Nmax = B->cmap.n;
1853: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C",
1854: "MatSeqDenseSetPreallocation_SeqDense",
1855: MatSeqDenseSetPreallocation_SeqDense);
1856: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C",
1857: "MatMatMult_SeqAIJ_SeqDense",
1858: MatMatMult_SeqAIJ_SeqDense);
1859: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",
1860: "MatMatMultSymbolic_SeqAIJ_SeqDense",
1861: MatMatMultSymbolic_SeqAIJ_SeqDense);
1862: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",
1863: "MatMatMultNumeric_SeqAIJ_SeqDense",
1864: MatMatMultNumeric_SeqAIJ_SeqDense);
1865: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
1866: return(0);
1867: }