Actual source code: dense.c
1: /*$Id: dense.c,v 1.208 2001/09/07 20:09:16 bsmith Exp $*/
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: #undef __FUNCT__
11: int MatAXPY_SeqDense(PetscScalar *alpha,Mat X,Mat Y,MatStructure str)
12: {
13: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
14: int N = X->m*X->n,one = 1;
17: BLaxpy_(&N,alpha,x->v,&one,y->v,&one);
18: PetscLogFlops(2*N-1);
19: return(0);
20: }
22: #undef __FUNCT__
24: int MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
25: {
26: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
27: int i,N = A->m*A->n,count = 0;
28: PetscScalar *v = a->v;
31: for (i=0; i<N; i++) {if (*v != 0.0) count++; v++;}
33: info->rows_global = (double)A->m;
34: info->columns_global = (double)A->n;
35: info->rows_local = (double)A->m;
36: info->columns_local = (double)A->n;
37: info->block_size = 1.0;
38: info->nz_allocated = (double)N;
39: info->nz_used = (double)count;
40: info->nz_unneeded = (double)(N-count);
41: info->assemblies = (double)A->num_ass;
42: info->mallocs = 0;
43: info->memory = A->mem;
44: info->fill_ratio_given = 0;
45: info->fill_ratio_needed = 0;
46: info->factor_mallocs = 0;
48: return(0);
49: }
51: #undef __FUNCT__
53: int MatScale_SeqDense(PetscScalar *alpha,Mat A)
54: {
55: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
56: int one = 1,nz;
59: nz = A->m*A->n;
60: BLscal_(&nz,alpha,a->v,&one);
61: PetscLogFlops(nz);
62: return(0);
63: }
64:
65: /* ---------------------------------------------------------------*/
66: /* COMMENT: I have chosen to hide row permutation in the pivots,
67: rather than put it in the Mat->row slot.*/
68: #undef __FUNCT__
70: int MatLUFactor_SeqDense(Mat A,IS row,IS col,MatLUInfo *minfo)
71: {
72: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
73: int info,ierr;
76: if (!mat->pivots) {
77: PetscMalloc((A->m+1)*sizeof(int),&mat->pivots);
78: PetscLogObjectMemory(A,A->m*sizeof(int));
79: }
80: A->factor = FACTOR_LU;
81: if (!A->m || !A->n) return(0);
82: #if defined(PETSC_MISSING_LAPACK_GETRF)
83: SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavilable.");
84: #else
85: LAgetrf_(&A->m,&A->n,mat->v,&A->m,mat->pivots,&info);
86: if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization");
87: if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
88: #endif
89: PetscLogFlops((2*A->n*A->n*A->n)/3);
90: return(0);
91: }
93: #undef __FUNCT__
95: int MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
96: {
97: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
98: int ierr;
99: Mat newi;
102: MatCreateSeqDense(A->comm,A->m,A->n,PETSC_NULL,&newi);
103: l = (Mat_SeqDense*)newi->data;
104: if (cpvalues == MAT_COPY_VALUES) {
105: PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));
106: }
107: newi->assembled = PETSC_TRUE;
108: *newmat = newi;
109: return(0);
110: }
112: #undef __FUNCT__
114: int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatLUInfo *info,Mat *fact)
115: {
119: MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);
120: return(0);
121: }
123: #undef __FUNCT__
125: int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact)
126: {
127: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l = (Mat_SeqDense*)(*fact)->data;
128: int ierr;
131: /* copy the numerical values */
132: PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));
133: (*fact)->factor = 0;
134: MatLUFactor(*fact,0,0,PETSC_NULL);
135: return(0);
136: }
138: #undef __FUNCT__
140: int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,PetscReal f,Mat *fact)
141: {
145: MatConvert(A,MATSAME,fact);
146: return(0);
147: }
149: #undef __FUNCT__
151: int MatCholeskyFactor_SeqDense(Mat A,IS perm,PetscReal f)
152: {
153: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
154: int info,ierr;
155:
157: if (mat->pivots) {
158: PetscFree(mat->pivots);
159: PetscLogObjectMemory(A,-A->m*sizeof(int));
160: mat->pivots = 0;
161: }
162: if (!A->m || !A->n) return(0);
163: #if defined(PETSC_MISSING_LAPACK_POTRF)
164: SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavilable.");
165: #else
166: LApotrf_("L",&A->n,mat->v,&A->m,&info);
167: if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %d",info-1);
168: #endif
169: A->factor = FACTOR_CHOLESKY;
170: PetscLogFlops((A->n*A->n*A->n)/3);
171: return(0);
172: }
174: #undef __FUNCT__
176: int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact)
177: {
181: MatCholeskyFactor_SeqDense(*fact,0,1.0);
182: return(0);
183: }
185: #undef __FUNCT__
187: int MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
188: {
189: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
190: int one = 1,info,ierr;
191: PetscScalar *x,*y;
192:
194: if (!A->m || !A->n) return(0);
195: VecGetArray(xx,&x);
196: VecGetArray(yy,&y);
197: PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
198: if (A->factor == FACTOR_LU) {
199: #if defined(PETSC_MISSING_LAPACK_GETRS)
200: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavilable.");
201: #else
202: LAgetrs_("N",&A->m,&one,mat->v,&A->m,mat->pivots,y,&A->m,&info);
203: if (info) SETERRQ(PETSC_ERR_LIB,"MBad solve");
204: #endif
205: } else if (A->factor == FACTOR_CHOLESKY){
206: #if defined(PETSC_MISSING_LAPACK_POTRS)
207: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavilable.");
208: #else
209: LApotrs_("L",&A->m,&one,mat->v,&A->m,y,&A->m,&info);
210: if (info) SETERRQ(PETSC_ERR_LIB,"MBad solve");
211: #endif
212: }
213: else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
214: VecRestoreArray(xx,&x);
215: VecRestoreArray(yy,&y);
216: PetscLogFlops(2*A->n*A->n - A->n);
217: return(0);
218: }
220: #undef __FUNCT__
222: int MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
223: {
224: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
225: int ierr,one = 1,info;
226: PetscScalar *x,*y;
227:
229: if (!A->m || !A->n) return(0);
230: VecGetArray(xx,&x);
231: VecGetArray(yy,&y);
232: PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
233: /* assume if pivots exist then use LU; else Cholesky */
234: if (mat->pivots) {
235: #if defined(PETSC_MISSING_LAPACK_GETRS)
236: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavilable.");
237: #else
238: LAgetrs_("T",&A->m,&one,mat->v,&A->m,mat->pivots,y,&A->m,&info);
239: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
240: #endif
241: } else {
242: #if defined(PETSC_MISSING_LAPACK_POTRS)
243: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavilable.");
244: #else
245: LApotrs_("L",&A->m,&one,mat->v,&A->m,y,&A->m,&info);
246: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
247: #endif
248: }
249: VecRestoreArray(xx,&x);
250: VecRestoreArray(yy,&y);
251: PetscLogFlops(2*A->n*A->n - A->n);
252: return(0);
253: }
255: #undef __FUNCT__
257: int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
258: {
259: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
260: int ierr,one = 1,info;
261: PetscScalar *x,*y,sone = 1.0;
262: Vec tmp = 0;
263:
265: VecGetArray(xx,&x);
266: VecGetArray(yy,&y);
267: if (!A->m || !A->n) return(0);
268: if (yy == zz) {
269: VecDuplicate(yy,&tmp);
270: PetscLogObjectParent(A,tmp);
271: VecCopy(yy,tmp);
272: }
273: PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
274: /* assume if pivots exist then use LU; else Cholesky */
275: if (mat->pivots) {
276: #if defined(PETSC_MISSING_LAPACK_GETRS)
277: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavilable.");
278: #else
279: LAgetrs_("N",&A->m,&one,mat->v,&A->m,mat->pivots,y,&A->m,&info);
280: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
281: #endif
282: } else {
283: #if defined(PETSC_MISSING_LAPACK_POTRS)
284: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavilable.");
285: #else
286: LApotrs_("L",&A->m,&one,mat->v,&A->m,y,&A->m,&info);
287: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
288: #endif
289: }
290: if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);}
291: else {VecAXPY(&sone,zz,yy);}
292: VecRestoreArray(xx,&x);
293: VecRestoreArray(yy,&y);
294: PetscLogFlops(2*A->n*A->n);
295: return(0);
296: }
298: #undef __FUNCT__
300: int MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
301: {
302: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
303: int one = 1,info,ierr;
304: PetscScalar *x,*y,sone = 1.0;
305: Vec tmp;
306:
308: if (!A->m || !A->n) return(0);
309: VecGetArray(xx,&x);
310: VecGetArray(yy,&y);
311: if (yy == zz) {
312: VecDuplicate(yy,&tmp);
313: PetscLogObjectParent(A,tmp);
314: VecCopy(yy,tmp);
315: }
316: PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
317: /* assume if pivots exist then use LU; else Cholesky */
318: if (mat->pivots) {
319: #if defined(PETSC_MISSING_LAPACK_GETRS)
320: SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavilable.");
321: #else
322: LAgetrs_("T",&A->m,&one,mat->v,&A->m,mat->pivots,y,&A->m,&info);
323: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
324: #endif
325: } else {
326: #if defined(PETSC_MISSING_LAPACK_POTRS)
327: SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavilable.");
328: #else
329: LApotrs_("L",&A->m,&one,mat->v,&A->m,y,&A->m,&info);
330: if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
331: #endif
332: }
333: if (tmp) {
334: VecAXPY(&sone,tmp,yy);
335: VecDestroy(tmp);
336: } else {
337: VecAXPY(&sone,zz,yy);
338: }
339: VecRestoreArray(xx,&x);
340: VecRestoreArray(yy,&y);
341: PetscLogFlops(2*A->n*A->n);
342: return(0);
343: }
344: /* ------------------------------------------------------------------*/
345: #undef __FUNCT__
347: int MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,
348: PetscReal shift,int its,int lits,Vec xx)
349: {
350: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
351: PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt;
352: int ierr,m = A->m,i;
353: #if !defined(PETSC_USE_COMPLEX)
354: int o = 1;
355: #endif
358: if (flag & SOR_ZERO_INITIAL_GUESS) {
359: /* this is a hack fix, should have another version without the second BLdot */
360: VecSet(&zero,xx);
361: }
362: VecGetArray(xx,&x);
363: VecGetArray(bb,&b);
364: its = its*lits;
365: if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits);
366: while (its--) {
367: if (flag & SOR_FORWARD_SWEEP){
368: for (i=0; i<m; i++) {
369: #if defined(PETSC_USE_COMPLEX)
370: /* cannot use BLAS dot for complex because compiler/linker is
371: not happy about returning a double complex */
372: int _i;
373: PetscScalar sum = b[i];
374: for (_i=0; _i<m; _i++) {
375: sum -= PetscConj(v[i+_i*m])*x[_i];
376: }
377: xt = sum;
378: #else
379: xt = b[i]-BLdot_(&m,v+i,&m,x,&o);
380: #endif
381: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
382: }
383: }
384: if (flag & SOR_BACKWARD_SWEEP) {
385: for (i=m-1; i>=0; i--) {
386: #if defined(PETSC_USE_COMPLEX)
387: /* cannot use BLAS dot for complex because compiler/linker is
388: not happy about returning a double complex */
389: int _i;
390: PetscScalar sum = b[i];
391: for (_i=0; _i<m; _i++) {
392: sum -= PetscConj(v[i+_i*m])*x[_i];
393: }
394: xt = sum;
395: #else
396: xt = b[i]-BLdot_(&m,v+i,&m,x,&o);
397: #endif
398: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
399: }
400: }
401: }
402: VecRestoreArray(bb,&b);
403: VecRestoreArray(xx,&x);
404: return(0);
405: }
407: /* -----------------------------------------------------------------*/
408: #undef __FUNCT__
410: int MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
411: {
412: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
413: PetscScalar *v = mat->v,*x,*y;
414: int ierr,_One=1;
415: PetscScalar _DOne=1.0,_DZero=0.0;
418: if (!A->m || !A->n) return(0);
419: VecGetArray(xx,&x);
420: VecGetArray(yy,&y);
421: LAgemv_("T",&(A->m),&(A->n),&_DOne,v,&(A->m),x,&_One,&_DZero,y,&_One);
422: VecRestoreArray(xx,&x);
423: VecRestoreArray(yy,&y);
424: PetscLogFlops(2*A->m*A->n - A->n);
425: return(0);
426: }
428: #undef __FUNCT__
430: int MatMult_SeqDense(Mat A,Vec xx,Vec yy)
431: {
432: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
433: PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
434: int ierr,_One=1;
437: if (!A->m || !A->n) return(0);
438: VecGetArray(xx,&x);
439: VecGetArray(yy,&y);
440: LAgemv_("N",&(A->m),&(A->n),&_DOne,v,&(A->m),x,&_One,&_DZero,y,&_One);
441: VecRestoreArray(xx,&x);
442: VecRestoreArray(yy,&y);
443: PetscLogFlops(2*A->m*A->n - A->m);
444: return(0);
445: }
447: #undef __FUNCT__
449: int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
450: {
451: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
452: PetscScalar *v = mat->v,*x,*y,_DOne=1.0;
453: int ierr,_One=1;
456: if (!A->m || !A->n) return(0);
457: if (zz != yy) {VecCopy(zz,yy);}
458: VecGetArray(xx,&x);
459: VecGetArray(yy,&y);
460: LAgemv_("N",&(A->m),&(A->n),&_DOne,v,&(A->m),x,&_One,&_DOne,y,&_One);
461: VecRestoreArray(xx,&x);
462: VecRestoreArray(yy,&y);
463: PetscLogFlops(2*A->m*A->n);
464: return(0);
465: }
467: #undef __FUNCT__
469: int MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
470: {
471: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
472: PetscScalar *v = mat->v,*x,*y;
473: int ierr,_One=1;
474: PetscScalar _DOne=1.0;
477: if (!A->m || !A->n) return(0);
478: if (zz != yy) {VecCopy(zz,yy);}
479: VecGetArray(xx,&x);
480: VecGetArray(yy,&y);
481: LAgemv_("T",&(A->m),&(A->n),&_DOne,v,&(A->m),x,&_One,&_DOne,y,&_One);
482: VecRestoreArray(xx,&x);
483: VecRestoreArray(yy,&y);
484: PetscLogFlops(2*A->m*A->n);
485: return(0);
486: }
488: /* -----------------------------------------------------------------*/
489: #undef __FUNCT__
491: int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals)
492: {
493: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
494: PetscScalar *v;
495: int i,ierr;
496:
498: *ncols = A->n;
499: if (cols) {
500: PetscMalloc((A->n+1)*sizeof(int),cols);
501: for (i=0; i<A->n; i++) (*cols)[i] = i;
502: }
503: if (vals) {
504: PetscMalloc((A->n+1)*sizeof(PetscScalar),vals);
505: v = mat->v + row;
506: for (i=0; i<A->n; i++) {(*vals)[i] = *v; v += A->m;}
507: }
508: return(0);
509: }
511: #undef __FUNCT__
513: int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals)
514: {
517: if (cols) {PetscFree(*cols);}
518: if (vals) {PetscFree(*vals); }
519: return(0);
520: }
521: /* ----------------------------------------------------------------*/
522: #undef __FUNCT__
524: int MatSetValues_SeqDense(Mat A,int m,int *indexm,int n,
525: int *indexn,PetscScalar *v,InsertMode addv)
526: {
527: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
528: int i,j;
529:
531: if (!mat->roworiented) {
532: if (addv == INSERT_VALUES) {
533: for (j=0; j<n; j++) {
534: if (indexn[j] < 0) {v += m; continue;}
535: #if defined(PETSC_USE_BOPT_g)
536: if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
537: #endif
538: for (i=0; i<m; i++) {
539: if (indexm[i] < 0) {v++; continue;}
540: #if defined(PETSC_USE_BOPT_g)
541: if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
542: #endif
543: mat->v[indexn[j]*A->m + indexm[i]] = *v++;
544: }
545: }
546: } else {
547: for (j=0; j<n; j++) {
548: if (indexn[j] < 0) {v += m; continue;}
549: #if defined(PETSC_USE_BOPT_g)
550: if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
551: #endif
552: for (i=0; i<m; i++) {
553: if (indexm[i] < 0) {v++; continue;}
554: #if defined(PETSC_USE_BOPT_g)
555: if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
556: #endif
557: mat->v[indexn[j]*A->m + indexm[i]] += *v++;
558: }
559: }
560: }
561: } else {
562: if (addv == INSERT_VALUES) {
563: for (i=0; i<m; i++) {
564: if (indexm[i] < 0) { v += n; continue;}
565: #if defined(PETSC_USE_BOPT_g)
566: if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
567: #endif
568: for (j=0; j<n; j++) {
569: if (indexn[j] < 0) { v++; continue;}
570: #if defined(PETSC_USE_BOPT_g)
571: if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
572: #endif
573: mat->v[indexn[j]*A->m + indexm[i]] = *v++;
574: }
575: }
576: } else {
577: for (i=0; i<m; i++) {
578: if (indexm[i] < 0) { v += n; continue;}
579: #if defined(PETSC_USE_BOPT_g)
580: if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
581: #endif
582: for (j=0; j<n; j++) {
583: if (indexn[j] < 0) { v++; continue;}
584: #if defined(PETSC_USE_BOPT_g)
585: if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
586: #endif
587: mat->v[indexn[j]*A->m + indexm[i]] += *v++;
588: }
589: }
590: }
591: }
592: return(0);
593: }
595: #undef __FUNCT__
597: int MatGetValues_SeqDense(Mat A,int m,int *indexm,int n,int *indexn,PetscScalar *v)
598: {
599: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
600: int i,j;
601: PetscScalar *vpt = v;
604: /* row-oriented output */
605: for (i=0; i<m; i++) {
606: for (j=0; j<n; j++) {
607: *vpt++ = mat->v[indexn[j]*A->m + indexm[i]];
608: }
609: }
610: return(0);
611: }
613: /* -----------------------------------------------------------------*/
615: #include petscsys.h
617: EXTERN_C_BEGIN
618: #undef __FUNCT__
620: int MatLoad_SeqDense(PetscViewer viewer,MatType type,Mat *A)
621: {
622: Mat_SeqDense *a;
623: Mat B;
624: int *scols,i,j,nz,ierr,fd,header[4],size;
625: int *rowlengths = 0,M,N,*cols;
626: PetscScalar *vals,*svals,*v,*w;
627: MPI_Comm comm = ((PetscObject)viewer)->comm;
630: MPI_Comm_size(comm,&size);
631: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor");
632: PetscViewerBinaryGetDescriptor(viewer,&fd);
633: PetscBinaryRead(fd,header,4,PETSC_INT);
634: if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
635: M = header[1]; N = header[2]; nz = header[3];
637: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
638: MatCreateSeqDense(comm,M,N,PETSC_NULL,A);
639: B = *A;
640: a = (Mat_SeqDense*)B->data;
641: v = a->v;
642: /* Allocate some temp space to read in the values and then flip them
643: from row major to column major */
644: PetscMalloc((M*N+1)*sizeof(PetscScalar),&w);
645: /* read in nonzero values */
646: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
647: /* now flip the values and store them in the matrix*/
648: for (j=0; j<N; j++) {
649: for (i=0; i<M; i++) {
650: *v++ =w[i*N+j];
651: }
652: }
653: PetscFree(w);
654: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
655: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
656: } else {
657: /* read row lengths */
658: PetscMalloc((M+1)*sizeof(int),&rowlengths);
659: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
661: /* create our matrix */
662: MatCreateSeqDense(comm,M,N,PETSC_NULL,A);
663: B = *A;
664: a = (Mat_SeqDense*)B->data;
665: v = a->v;
667: /* read column indices and nonzeros */
668: PetscMalloc((nz+1)*sizeof(int),&scols);
669: cols = scols;
670: PetscBinaryRead(fd,cols,nz,PETSC_INT);
671: PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);
672: vals = svals;
673: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
675: /* insert into matrix */
676: for (i=0; i<M; i++) {
677: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
678: svals += rowlengths[i]; scols += rowlengths[i];
679: }
680: PetscFree(vals);
681: PetscFree(cols);
682: PetscFree(rowlengths);
684: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
685: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
686: }
687: return(0);
688: }
689: EXTERN_C_END
691: #include petscsys.h
693: #undef __FUNCT__
695: static int MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
696: {
697: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
698: int ierr,i,j;
699: char *name;
700: PetscScalar *v;
701: PetscViewerFormat format;
704: PetscObjectGetName((PetscObject)A,&name);
705: PetscViewerGetFormat(viewer,&format);
706: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
707: return(0); /* do nothing for now */
708: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
709: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
710: for (i=0; i<A->m; i++) {
711: v = a->v + i;
712: PetscViewerASCIIPrintf(viewer,"row %d:",i);
713: for (j=0; j<A->n; j++) {
714: #if defined(PETSC_USE_COMPLEX)
715: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
716: PetscViewerASCIIPrintf(viewer," (%d, %g + %g i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));
717: } else if (PetscRealPart(*v)) {
718: PetscViewerASCIIPrintf(viewer," (%d, %g) ",j,PetscRealPart(*v));
719: }
720: #else
721: if (*v) {
722: PetscViewerASCIIPrintf(viewer," (%d, %g) ",j,*v);
723: }
724: #endif
725: v += A->m;
726: }
727: PetscViewerASCIIPrintf(viewer,"n");
728: }
729: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
730: } else {
731: PetscViewerASCIIUseTabs(viewer,PETSC_NO);
732: #if defined(PETSC_USE_COMPLEX)
733: PetscTruth allreal = PETSC_TRUE;
734: /* determine if matrix has all real values */
735: v = a->v;
736: for (i=0; i<A->m*A->n; i++) {
737: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;}
738: }
739: #endif
740: if (format == PETSC_VIEWER_ASCII_MATLAB) {
741: PetscObjectGetName((PetscObject)A,&name);
742: PetscViewerASCIIPrintf(viewer,"%% Size = %d %d n",A->m,A->n);
743: PetscViewerASCIIPrintf(viewer,"%s = zeros(%d,%d);n",name,A->m,A->n);
744: PetscViewerASCIIPrintf(viewer,"%s = [n",name);
745: }
747: for (i=0; i<A->m; i++) {
748: v = a->v + i;
749: for (j=0; j<A->n; j++) {
750: #if defined(PETSC_USE_COMPLEX)
751: if (allreal) {
752: PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v)); v += A->m;
753: } else {
754: PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v)); v += A->m;
755: }
756: #else
757: PetscViewerASCIIPrintf(viewer,"%6.4e ",*v); v += A->m;
758: #endif
759: }
760: PetscViewerASCIIPrintf(viewer,"n");
761: }
762: if (format == PETSC_VIEWER_ASCII_MATLAB) {
763: PetscViewerASCIIPrintf(viewer,"];n");
764: }
765: PetscViewerASCIIUseTabs(viewer,PETSC_YES);
766: }
767: PetscViewerFlush(viewer);
768: return(0);
769: }
771: #undef __FUNCT__
773: static int MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
774: {
775: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
776: int ict,j,n = A->n,m = A->m,i,fd,*col_lens,ierr,nz = m*n;
777: PetscScalar *v,*anonz,*vals;
778: PetscViewerFormat format;
779:
781: PetscViewerBinaryGetDescriptor(viewer,&fd);
783: PetscViewerGetFormat(viewer,&format);
784: if (format == PETSC_VIEWER_BINARY_NATIVE) {
785: /* store the matrix as a dense matrix */
786: PetscMalloc(4*sizeof(int),&col_lens);
787: col_lens[0] = MAT_FILE_COOKIE;
788: col_lens[1] = m;
789: col_lens[2] = n;
790: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
791: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);
792: PetscFree(col_lens);
794: /* write out matrix, by rows */
795: PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);
796: v = a->v;
797: for (i=0; i<m; i++) {
798: for (j=0; j<n; j++) {
799: vals[i + j*m] = *v++;
800: }
801: }
802: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);
803: PetscFree(vals);
804: } else {
805: PetscMalloc((4+nz)*sizeof(int),&col_lens);
806: col_lens[0] = MAT_FILE_COOKIE;
807: col_lens[1] = m;
808: col_lens[2] = n;
809: col_lens[3] = nz;
811: /* store lengths of each row and write (including header) to file */
812: for (i=0; i<m; i++) col_lens[4+i] = n;
813: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);
815: /* Possibly should write in smaller increments, not whole matrix at once? */
816: /* store column indices (zero start index) */
817: ict = 0;
818: for (i=0; i<m; i++) {
819: for (j=0; j<n; j++) col_lens[ict++] = j;
820: }
821: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);
822: PetscFree(col_lens);
824: /* store nonzero values */
825: PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);
826: ict = 0;
827: for (i=0; i<m; i++) {
828: v = a->v + i;
829: for (j=0; j<n; j++) {
830: anonz[ict++] = *v; v += A->m;
831: }
832: }
833: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);
834: PetscFree(anonz);
835: }
836: return(0);
837: }
839: #undef __FUNCT__
841: int MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
842: {
843: Mat A = (Mat) Aa;
844: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
845: int m = A->m,n = A->n,color,i,j,ierr;
846: PetscScalar *v = a->v;
847: PetscViewer viewer;
848: PetscDraw popup;
849: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
850: PetscViewerFormat format;
854: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
855: PetscViewerGetFormat(viewer,&format);
856: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
858: /* Loop over matrix elements drawing boxes */
859: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
860: /* Blue for negative and Red for positive */
861: color = PETSC_DRAW_BLUE;
862: for(j = 0; j < n; j++) {
863: x_l = j;
864: x_r = x_l + 1.0;
865: for(i = 0; i < m; i++) {
866: y_l = m - i - 1.0;
867: y_r = y_l + 1.0;
868: #if defined(PETSC_USE_COMPLEX)
869: if (PetscRealPart(v[j*m+i]) > 0.) {
870: color = PETSC_DRAW_RED;
871: } else if (PetscRealPart(v[j*m+i]) < 0.) {
872: color = PETSC_DRAW_BLUE;
873: } else {
874: continue;
875: }
876: #else
877: if (v[j*m+i] > 0.) {
878: color = PETSC_DRAW_RED;
879: } else if (v[j*m+i] < 0.) {
880: color = PETSC_DRAW_BLUE;
881: } else {
882: continue;
883: }
884: #endif
885: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
886: }
887: }
888: } else {
889: /* use contour shading to indicate magnitude of values */
890: /* first determine max of all nonzero values */
891: for(i = 0; i < m*n; i++) {
892: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
893: }
894: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
895: ierr = PetscDrawGetPopup(draw,&popup);
896: if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);}
897: for(j = 0; j < n; j++) {
898: x_l = j;
899: x_r = x_l + 1.0;
900: for(i = 0; i < m; i++) {
901: y_l = m - i - 1.0;
902: y_r = y_l + 1.0;
903: color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
904: ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
905: }
906: }
907: }
908: return(0);
909: }
911: #undef __FUNCT__
913: int MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
914: {
915: PetscDraw draw;
916: PetscTruth isnull;
917: PetscReal xr,yr,xl,yl,h,w;
918: int ierr;
921: PetscViewerDrawGetDraw(viewer,0,&draw);
922: PetscDrawIsNull(draw,&isnull);
923: if (isnull == PETSC_TRUE) return(0);
925: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
926: xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0;
927: xr += w; yr += h; xl = -w; yl = -h;
928: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
929: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
930: PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
931: return(0);
932: }
934: #undef __FUNCT__
936: int MatView_SeqDense(Mat A,PetscViewer viewer)
937: {
938: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
939: int ierr;
940: PetscTruth issocket,isascii,isbinary,isdraw;
943: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);
944: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
945: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
946: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
948: if (issocket) {
949: PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);
950: } else if (isascii) {
951: MatView_SeqDense_ASCII(A,viewer);
952: } else if (isbinary) {
953: MatView_SeqDense_Binary(A,viewer);
954: } else if (isdraw) {
955: MatView_SeqDense_Draw(A,viewer);
956: } else {
957: SETERRQ1(1,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name);
958: }
959: return(0);
960: }
962: #undef __FUNCT__
964: int MatDestroy_SeqDense(Mat mat)
965: {
966: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
967: int ierr;
970: #if defined(PETSC_USE_LOG)
971: PetscLogObjectState((PetscObject)mat,"Rows %d Cols %d",mat->m,mat->n);
972: #endif
973: if (l->pivots) {PetscFree(l->pivots);}
974: if (!l->user_alloc) {PetscFree(l->v);}
975: PetscFree(l);
976: return(0);
977: }
979: #undef __FUNCT__
981: int MatTranspose_SeqDense(Mat A,Mat *matout)
982: {
983: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
984: int k,j,m,n,ierr;
985: PetscScalar *v,tmp;
988: v = mat->v; m = A->m; n = A->n;
989: if (!matout) { /* in place transpose */
990: if (m != n) { /* malloc temp to hold transpose */
991: PetscScalar *w;
992: PetscMalloc((m+1)*(n+1)*sizeof(PetscScalar),&w);
993: for (j=0; j<m; j++) {
994: for (k=0; k<n; k++) {
995: w[k + j*n] = v[j + k*m];
996: }
997: }
998: PetscMemcpy(v,w,m*n*sizeof(PetscScalar));
999: PetscFree(w);
1000: } else {
1001: for (j=0; j<m; j++) {
1002: for (k=0; k<j; k++) {
1003: tmp = v[j + k*n];
1004: v[j + k*n] = v[k + j*n];
1005: v[k + j*n] = tmp;
1006: }
1007: }
1008: }
1009: } else { /* out-of-place transpose */
1010: Mat tmat;
1011: Mat_SeqDense *tmatd;
1012: PetscScalar *v2;
1014: ierr = MatCreateSeqDense(A->comm,A->n,A->m,PETSC_NULL,&tmat);
1015: tmatd = (Mat_SeqDense*)tmat->data;
1016: v = mat->v; v2 = tmatd->v;
1017: for (j=0; j<n; j++) {
1018: for (k=0; k<m; k++) v2[j + k*n] = v[k + j*m];
1019: }
1020: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1021: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1022: *matout = tmat;
1023: }
1024: return(0);
1025: }
1027: #undef __FUNCT__
1029: int MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg)
1030: {
1031: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1032: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1033: int ierr,i;
1034: PetscScalar *v1 = mat1->v,*v2 = mat2->v;
1035: PetscTruth flag;
1038: PetscTypeCompare((PetscObject)A2,MATSEQDENSE,&flag);
1039: if (!flag) SETERRQ(PETSC_ERR_SUP,"Matrices should be of type MATSEQDENSE");
1040: if (A1->m != A2->m) {*flg = PETSC_FALSE; return(0);}
1041: if (A1->n != A2->n) {*flg = PETSC_FALSE; return(0);}
1042: for (i=0; i<A1->m*A1->n; i++) {
1043: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1044: v1++; v2++;
1045: }
1046: *flg = PETSC_TRUE;
1047: return(0);
1048: }
1050: #undef __FUNCT__
1052: int MatGetDiagonal_SeqDense(Mat A,Vec v)
1053: {
1054: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1055: int ierr,i,n,len;
1056: PetscScalar *x,zero = 0.0;
1059: VecSet(&zero,v);
1060: VecGetSize(v,&n);
1061: VecGetArray(v,&x);
1062: len = PetscMin(A->m,A->n);
1063: if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1064: for (i=0; i<len; i++) {
1065: x[i] = mat->v[i*A->m + i];
1066: }
1067: VecRestoreArray(v,&x);
1068: return(0);
1069: }
1071: #undef __FUNCT__
1073: int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1074: {
1075: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1076: PetscScalar *l,*r,x,*v;
1077: int ierr,i,j,m = A->m,n = A->n;
1080: if (ll) {
1081: VecGetSize(ll,&m);
1082: VecGetArray(ll,&l);
1083: if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1084: for (i=0; i<m; i++) {
1085: x = l[i];
1086: v = mat->v + i;
1087: for (j=0; j<n; j++) { (*v) *= x; v+= m;}
1088: }
1089: VecRestoreArray(ll,&l);
1090: PetscLogFlops(n*m);
1091: }
1092: if (rr) {
1093: VecGetSize(rr,&n);
1094: VecGetArray(rr,&r);
1095: if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1096: for (i=0; i<n; i++) {
1097: x = r[i];
1098: v = mat->v + i*m;
1099: for (j=0; j<m; j++) { (*v++) *= x;}
1100: }
1101: VecRestoreArray(rr,&r);
1102: PetscLogFlops(n*m);
1103: }
1104: return(0);
1105: }
1107: #undef __FUNCT__
1109: int MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1110: {
1111: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1112: PetscScalar *v = mat->v;
1113: PetscReal sum = 0.0;
1114: int i,j;
1117: if (type == NORM_FROBENIUS) {
1118: for (i=0; i<A->n*A->m; i++) {
1119: #if defined(PETSC_USE_COMPLEX)
1120: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1121: #else
1122: sum += (*v)*(*v); v++;
1123: #endif
1124: }
1125: *nrm = sqrt(sum);
1126: PetscLogFlops(2*A->n*A->m);
1127: } else if (type == NORM_1) {
1128: *nrm = 0.0;
1129: for (j=0; j<A->n; j++) {
1130: sum = 0.0;
1131: for (i=0; i<A->m; i++) {
1132: sum += PetscAbsScalar(*v); v++;
1133: }
1134: if (sum > *nrm) *nrm = sum;
1135: }
1136: PetscLogFlops(A->n*A->m);
1137: } else if (type == NORM_INFINITY) {
1138: *nrm = 0.0;
1139: for (j=0; j<A->m; j++) {
1140: v = mat->v + j;
1141: sum = 0.0;
1142: for (i=0; i<A->n; i++) {
1143: sum += PetscAbsScalar(*v); v += A->m;
1144: }
1145: if (sum > *nrm) *nrm = sum;
1146: }
1147: PetscLogFlops(A->n*A->m);
1148: } else {
1149: SETERRQ(PETSC_ERR_SUP,"No two norm");
1150: }
1151: return(0);
1152: }
1154: #undef __FUNCT__
1156: int MatSetOption_SeqDense(Mat A,MatOption op)
1157: {
1158: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1159:
1161: switch (op) {
1162: case MAT_ROW_ORIENTED:
1163: aij->roworiented = PETSC_TRUE;
1164: break;
1165: case MAT_COLUMN_ORIENTED:
1166: aij->roworiented = PETSC_FALSE;
1167: break;
1168: case MAT_ROWS_SORTED:
1169: case MAT_ROWS_UNSORTED:
1170: case MAT_COLUMNS_SORTED:
1171: case MAT_COLUMNS_UNSORTED:
1172: case MAT_NO_NEW_NONZERO_LOCATIONS:
1173: case MAT_YES_NEW_NONZERO_LOCATIONS:
1174: case MAT_NEW_NONZERO_LOCATION_ERR:
1175: case MAT_NO_NEW_DIAGONALS:
1176: case MAT_YES_NEW_DIAGONALS:
1177: case MAT_IGNORE_OFF_PROC_ENTRIES:
1178: case MAT_USE_HASH_TABLE:
1179: PetscLogInfo(A,"MatSetOption_SeqDense:Option ignoredn");
1180: break;
1181: default:
1182: SETERRQ(PETSC_ERR_SUP,"unknown option");
1183: }
1184: return(0);
1185: }
1187: #undef __FUNCT__
1189: int MatZeroEntries_SeqDense(Mat A)
1190: {
1191: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1192: int ierr;
1195: PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));
1196: return(0);
1197: }
1199: #undef __FUNCT__
1201: int MatGetBlockSize_SeqDense(Mat A,int *bs)
1202: {
1204: *bs = 1;
1205: return(0);
1206: }
1208: #undef __FUNCT__
1210: int MatZeroRows_SeqDense(Mat A,IS is,PetscScalar *diag)
1211: {
1212: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1213: int n = A->n,i,j,ierr,N,*rows;
1214: PetscScalar *slot;
1217: ISGetLocalSize(is,&N);
1218: ISGetIndices(is,&rows);
1219: for (i=0; i<N; i++) {
1220: slot = l->v + rows[i];
1221: for (j=0; j<n; j++) { *slot = 0.0; slot += n;}
1222: }
1223: if (diag) {
1224: for (i=0; i<N; i++) {
1225: slot = l->v + (n+1)*rows[i];
1226: *slot = *diag;
1227: }
1228: }
1229: ISRestoreIndices(is,&rows);
1230: return(0);
1231: }
1233: #undef __FUNCT__
1235: int MatGetArray_SeqDense(Mat A,PetscScalar **array)
1236: {
1237: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1240: *array = mat->v;
1241: return(0);
1242: }
1244: #undef __FUNCT__
1246: int MatRestoreArray_SeqDense(Mat A,PetscScalar **array)
1247: {
1249: *array = 0; /* user cannot accidently use the array later */
1250: return(0);
1251: }
1253: #undef __FUNCT__
1255: static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B)
1256: {
1257: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1258: int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols;
1259: PetscScalar *av,*bv,*v = mat->v;
1260: Mat newmat;
1263: ISGetIndices(isrow,&irow);
1264: ISGetIndices(iscol,&icol);
1265: ISGetLocalSize(isrow,&nrows);
1266: ISGetLocalSize(iscol,&ncols);
1267:
1268: /* Check submatrixcall */
1269: if (scall == MAT_REUSE_MATRIX) {
1270: int n_cols,n_rows;
1271: MatGetSize(*B,&n_rows,&n_cols);
1272: if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
1273: newmat = *B;
1274: } else {
1275: /* Create and fill new matrix */
1276: MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);
1277: }
1279: /* Now extract the data pointers and do the copy,column at a time */
1280: bv = ((Mat_SeqDense*)newmat->data)->v;
1281:
1282: for (i=0; i<ncols; i++) {
1283: av = v + m*icol[i];
1284: for (j=0; j<nrows; j++) {
1285: *bv++ = av[irow[j]];
1286: }
1287: }
1289: /* Assemble the matrices so that the correct flags are set */
1290: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1291: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1293: /* Free work space */
1294: ISRestoreIndices(isrow,&irow);
1295: ISRestoreIndices(iscol,&icol);
1296: *B = newmat;
1297: return(0);
1298: }
1300: #undef __FUNCT__
1302: int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B)
1303: {
1304: int ierr,i;
1307: if (scall == MAT_INITIAL_MATRIX) {
1308: PetscMalloc((n+1)*sizeof(Mat),B);
1309: }
1311: for (i=0; i<n; i++) {
1312: MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1313: }
1314: return(0);
1315: }
1317: #undef __FUNCT__
1319: int MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1320: {
1321: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data;
1322: int ierr;
1323: PetscTruth flag;
1326: PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flag);
1327: if (!flag) {
1328: MatCopy_Basic(A,B,str);
1329: return(0);
1330: }
1331: if (A->m != B->m || A->n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1332: PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));
1333: return(0);
1334: }
1336: #undef __FUNCT__
1338: int MatSetUpPreallocation_SeqDense(Mat A)
1339: {
1340: int ierr;
1343: MatSeqDenseSetPreallocation(A,0);
1344: return(0);
1345: }
1347: /* -------------------------------------------------------------------*/
1348: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
1349: MatGetRow_SeqDense,
1350: MatRestoreRow_SeqDense,
1351: MatMult_SeqDense,
1352: MatMultAdd_SeqDense,
1353: MatMultTranspose_SeqDense,
1354: MatMultTransposeAdd_SeqDense,
1355: MatSolve_SeqDense,
1356: MatSolveAdd_SeqDense,
1357: MatSolveTranspose_SeqDense,
1358: MatSolveTransposeAdd_SeqDense,
1359: MatLUFactor_SeqDense,
1360: MatCholeskyFactor_SeqDense,
1361: MatRelax_SeqDense,
1362: MatTranspose_SeqDense,
1363: MatGetInfo_SeqDense,
1364: MatEqual_SeqDense,
1365: MatGetDiagonal_SeqDense,
1366: MatDiagonalScale_SeqDense,
1367: MatNorm_SeqDense,
1368: 0,
1369: 0,
1370: 0,
1371: MatSetOption_SeqDense,
1372: MatZeroEntries_SeqDense,
1373: MatZeroRows_SeqDense,
1374: MatLUFactorSymbolic_SeqDense,
1375: MatLUFactorNumeric_SeqDense,
1376: MatCholeskyFactorSymbolic_SeqDense,
1377: MatCholeskyFactorNumeric_SeqDense,
1378: MatSetUpPreallocation_SeqDense,
1379: 0,
1380: 0,
1381: MatGetArray_SeqDense,
1382: MatRestoreArray_SeqDense,
1383: MatDuplicate_SeqDense,
1384: 0,
1385: 0,
1386: 0,
1387: 0,
1388: MatAXPY_SeqDense,
1389: MatGetSubMatrices_SeqDense,
1390: 0,
1391: MatGetValues_SeqDense,
1392: MatCopy_SeqDense,
1393: 0,
1394: MatScale_SeqDense,
1395: 0,
1396: 0,
1397: 0,
1398: MatGetBlockSize_SeqDense,
1399: 0,
1400: 0,
1401: 0,
1402: 0,
1403: 0,
1404: 0,
1405: 0,
1406: 0,
1407: 0,
1408: 0,
1409: MatDestroy_SeqDense,
1410: MatView_SeqDense,
1411: MatGetPetscMaps_Petsc};
1413: #undef __FUNCT__
1415: /*@C
1416: MatCreateSeqDense - Creates a sequential dense matrix that
1417: is stored in column major order (the usual Fortran 77 manner). Many
1418: of the matrix operations use the BLAS and LAPACK routines.
1420: Collective on MPI_Comm
1422: Input Parameters:
1423: + comm - MPI communicator, set to PETSC_COMM_SELF
1424: . m - number of rows
1425: . n - number of columns
1426: - data - optional location of matrix data. Set data=PETSC_NULL for PETSc
1427: to control all matrix memory allocation.
1429: Output Parameter:
1430: . A - the matrix
1432: Notes:
1433: The data input variable is intended primarily for Fortran programmers
1434: who wish to allocate their own matrix memory space. Most users should
1435: set data=PETSC_NULL.
1437: Level: intermediate
1439: .keywords: dense, matrix, LAPACK, BLAS
1441: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
1442: @*/
1443: int MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A)
1444: {
1448: MatCreate(comm,m,n,m,n,A);
1449: MatSetType(*A,MATSEQDENSE);
1450: MatSeqDenseSetPreallocation(*A,data);
1451: return(0);
1452: }
1454: #undef __FUNCT__
1456: /*@C
1457: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
1459: Collective on MPI_Comm
1461: Input Parameters:
1462: + A - the matrix
1463: - data - the array (or PETSC_NULL)
1465: Notes:
1466: The data input variable is intended primarily for Fortran programmers
1467: who wish to allocate their own matrix memory space. Most users should
1468: set data=PETSC_NULL.
1470: Level: intermediate
1472: .keywords: dense, matrix, LAPACK, BLAS
1474: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
1475: @*/
1476: int MatSeqDenseSetPreallocation(Mat B,PetscScalar *data)
1477: {
1478: Mat_SeqDense *b;
1479: int ierr;
1480: PetscTruth flg2;
1483: PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flg2);
1484: if (!flg2) return(0);
1485: B->preallocated = PETSC_TRUE;
1486: b = (Mat_SeqDense*)B->data;
1487: if (!data) {
1488: ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);
1489: ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));
1490: b->user_alloc = PETSC_FALSE;
1491: PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar));
1492: } else { /* user-allocated storage */
1493: b->v = data;
1494: b->user_alloc = PETSC_TRUE;
1495: }
1496: return(0);
1497: }
1499: EXTERN_C_BEGIN
1500: #undef __FUNCT__
1502: int MatCreate_SeqDense(Mat B)
1503: {
1504: Mat_SeqDense *b;
1505: int ierr,size;
1508: MPI_Comm_size(B->comm,&size);
1509: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1511: B->m = B->M = PetscMax(B->m,B->M);
1512: B->n = B->N = PetscMax(B->n,B->N);
1514: ierr = PetscNew(Mat_SeqDense,&b);
1515: ierr = PetscMemzero(b,sizeof(Mat_SeqDense));
1516: ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
1517: B->factor = 0;
1518: B->mapping = 0;
1519: PetscLogObjectMemory(B,sizeof(struct _p_Mat));
1520: B->data = (void*)b;
1522: PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);
1523: PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);
1525: b->pivots = 0;
1526: b->roworiented = PETSC_TRUE;
1527: b->v = 0;
1529: return(0);
1530: }
1531: EXTERN_C_END