Actual source code: bdiag.c
1: /*$Id: bdiag.c,v 1.198 2001/08/07 03:02:53 balay Exp $*/
3: /* Block diagonal matrix format */
5: #include src/mat/impls/bdiag/seq/bdiag.h
6: #include src/vec/vecimpl.h
7: #include src/inline/ilu.h
9: EXTERN int MatSetValues_SeqBDiag_1(Mat,int,int *,int,int *,PetscScalar *,InsertMode);
10: EXTERN int MatSetValues_SeqBDiag_N(Mat,int,int *,int,int *,PetscScalar *,InsertMode);
11: EXTERN int MatGetValues_SeqBDiag_1(Mat,int,int *,int,int *,PetscScalar *);
12: EXTERN int MatGetValues_SeqBDiag_N(Mat,int,int *,int,int *,PetscScalar *);
13: EXTERN int MatMult_SeqBDiag_1(Mat,Vec,Vec);
14: EXTERN int MatMult_SeqBDiag_2(Mat,Vec,Vec);
15: EXTERN int MatMult_SeqBDiag_3(Mat,Vec,Vec);
16: EXTERN int MatMult_SeqBDiag_4(Mat,Vec,Vec);
17: EXTERN int MatMult_SeqBDiag_5(Mat,Vec,Vec);
18: EXTERN int MatMult_SeqBDiag_N(Mat,Vec,Vec);
19: EXTERN int MatMultAdd_SeqBDiag_1(Mat,Vec,Vec,Vec);
20: EXTERN int MatMultAdd_SeqBDiag_2(Mat,Vec,Vec,Vec);
21: EXTERN int MatMultAdd_SeqBDiag_3(Mat,Vec,Vec,Vec);
22: EXTERN int MatMultAdd_SeqBDiag_4(Mat,Vec,Vec,Vec);
23: EXTERN int MatMultAdd_SeqBDiag_5(Mat,Vec,Vec,Vec);
24: EXTERN int MatMultAdd_SeqBDiag_N(Mat,Vec,Vec,Vec);
25: EXTERN int MatMultTranspose_SeqBDiag_1(Mat,Vec,Vec);
26: EXTERN int MatMultTranspose_SeqBDiag_N(Mat,Vec,Vec);
27: EXTERN int MatMultTransposeAdd_SeqBDiag_1(Mat,Vec,Vec,Vec);
28: EXTERN int MatMultTransposeAdd_SeqBDiag_N(Mat,Vec,Vec,Vec);
29: EXTERN int MatRelax_SeqBDiag_N(Mat,Vec,PetscReal,MatSORType,PetscReal,int,int,Vec);
30: EXTERN int MatRelax_SeqBDiag_1(Mat,Vec,PetscReal,MatSORType,PetscReal,int,int,Vec);
31: EXTERN int MatView_SeqBDiag(Mat,PetscViewer);
32: EXTERN int MatGetInfo_SeqBDiag(Mat,MatInfoType,MatInfo*);
33: EXTERN int MatGetRow_SeqBDiag(Mat,int,int *,int **,PetscScalar **);
34: EXTERN int MatRestoreRow_SeqBDiag(Mat,int,int *,int **,PetscScalar **);
35: EXTERN int MatTranspose_SeqBDiag(Mat,Mat *);
36: EXTERN int MatNorm_SeqBDiag(Mat,NormType,PetscReal *);
38: #undef __FUNCT__
40: int MatDestroy_SeqBDiag(Mat A)
41: {
42: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
43: int i,bs = a->bs,ierr;
46: #if defined(PETSC_USE_LOG)
47: PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d, BSize=%d, NDiag=%d",A->m,A->n,a->nz,a->bs,a->nd);
48: #endif
49: if (!a->user_alloc) { /* Free the actual diagonals */
50: for (i=0; i<a->nd; i++) {
51: if (a->diag[i] > 0) {
52: PetscFree(a->diagv[i] + bs*bs*a->diag[i]);
53: } else {
54: PetscFree(a->diagv[i]);
55: }
56: }
57: }
58: if (a->pivot) {PetscFree(a->pivot);}
59: PetscFree(a->diagv);
60: PetscFree(a->diag);
61: PetscFree(a->colloc);
62: PetscFree(a->dvalue);
63: PetscFree(a);
64: return(0);
65: }
67: #undef __FUNCT__
69: int MatAssemblyEnd_SeqBDiag(Mat A,MatAssemblyType mode)
70: {
71: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
72: int i,k,temp,*diag = a->diag,*bdlen = a->bdlen;
73: PetscScalar *dtemp,**dv = a->diagv;
76: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
78: /* Sort diagonals */
79: for (i=0; i<a->nd; i++) {
80: for (k=i+1; k<a->nd; k++) {
81: if (diag[i] < diag[k]) {
82: temp = diag[i];
83: diag[i] = diag[k];
84: diag[k] = temp;
85: temp = bdlen[i];
86: bdlen[i] = bdlen[k];
87: bdlen[k] = temp;
88: dtemp = dv[i];
89: dv[i] = dv[k];
90: dv[k] = dtemp;
91: }
92: }
93: }
95: /* Set location of main diagonal */
96: for (i=0; i<a->nd; i++) {
97: if (!a->diag[i]) {a->mainbd = i; break;}
98: }
99: PetscLogInfo(A,"MatAssemblyEnd_SeqBDiag:Number diagonals %d,memory used %d, block size %dn",a->nd,a->maxnz,a->bs);
100: return(0);
101: }
103: #undef __FUNCT__
105: int MatSetOption_SeqBDiag(Mat A,MatOption op)
106: {
107: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
110: switch (op) {
111: case MAT_NO_NEW_NONZERO_LOCATIONS:
112: a->nonew = 1;
113: break;
114: case MAT_YES_NEW_NONZERO_LOCATIONS:
115: a->nonew = 0;
116: break;
117: case MAT_NO_NEW_DIAGONALS:
118: a->nonew_diag = 1;
119: break;
120: case MAT_YES_NEW_DIAGONALS:
121: a->nonew_diag = 0;
122: break;
123: case MAT_COLUMN_ORIENTED:
124: a->roworiented = PETSC_FALSE;
125: break;
126: case MAT_ROW_ORIENTED:
127: a->roworiented = PETSC_TRUE;
128: break;
129: case MAT_ROWS_SORTED:
130: case MAT_ROWS_UNSORTED:
131: case MAT_COLUMNS_SORTED:
132: case MAT_COLUMNS_UNSORTED:
133: case MAT_IGNORE_OFF_PROC_ENTRIES:
134: case MAT_NEW_NONZERO_LOCATION_ERR:
135: case MAT_NEW_NONZERO_ALLOCATION_ERR:
136: case MAT_USE_HASH_TABLE:
137: PetscLogInfo(A,"MatSetOption_SeqBDiag:Option ignoredn");
138: break;
139: default:
140: SETERRQ(PETSC_ERR_SUP,"unknown option");
141: }
142: return(0);
143: }
145: #undef __FUNCT__
147: int MatPrintHelp_SeqBDiag(Mat A)
148: {
149: static PetscTruth called = PETSC_FALSE;
150: MPI_Comm comm = A->comm;
151: int ierr;
154: if (called) {return(0);} else called = PETSC_TRUE;
155: (*PetscHelpPrintf)(comm," Options for MATSEQBDIAG and MATMPIBDIAG matrix formats:n");
156: (*PetscHelpPrintf)(comm," -mat_block_size <block_size>n");
157: (*PetscHelpPrintf)(comm," -mat_bdiag_diags <d1,d2,d3,...> (diagonal numbers)n");
158: (*PetscHelpPrintf)(comm," (for example) -mat_bdiag_diags -5,-1,0,1,5n");
159: return(0);
160: }
162: #undef __FUNCT__
164: static int MatGetDiagonal_SeqBDiag_N(Mat A,Vec v)
165: {
166: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
167: int ierr,i,j,n,len,ibase,bs = a->bs,iloc;
168: PetscScalar *x,*dd,zero = 0.0;
171: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
172: VecSet(&zero,v);
173: VecGetLocalSize(v,&n);
174: if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
175: if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal not set");
176: len = PetscMin(a->mblock,a->nblock);
177: dd = a->diagv[a->mainbd];
178: VecGetArray(v,&x);
179: for (i=0; i<len; i++) {
180: ibase = i*bs*bs; iloc = i*bs;
181: for (j=0; j<bs; j++) x[j + iloc] = dd[ibase + j*(bs+1)];
182: }
183: VecRestoreArray(v,&x);
184: return(0);
185: }
187: #undef __FUNCT__
189: static int MatGetDiagonal_SeqBDiag_1(Mat A,Vec v)
190: {
191: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
192: int ierr,i,n,len;
193: PetscScalar *x,*dd,zero = 0.0;
196: VecSet(&zero,v);
197: VecGetLocalSize(v,&n);
198: if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
199: if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal not set");
200: dd = a->diagv[a->mainbd];
201: len = PetscMin(A->m,A->n);
202: VecGetArray(v,&x);
203: for (i=0; i<len; i++) x[i] = dd[i];
204: VecRestoreArray(v,&x);
205: return(0);
206: }
208: #undef __FUNCT__
210: int MatZeroEntries_SeqBDiag(Mat A)
211: {
212: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
213: int d,i,len,bs = a->bs;
214: PetscScalar *dv;
217: for (d=0; d<a->nd; d++) {
218: dv = a->diagv[d];
219: if (a->diag[d] > 0) {
220: dv += bs*bs*a->diag[d];
221: }
222: len = a->bdlen[d]*bs*bs;
223: for (i=0; i<len; i++) dv[i] = 0.0;
224: }
225: return(0);
226: }
228: #undef __FUNCT__
230: int MatGetBlockSize_SeqBDiag(Mat A,int *bs)
231: {
232: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
235: *bs = a->bs;
236: return(0);
237: }
239: #undef __FUNCT__
241: int MatZeroRows_SeqBDiag(Mat A,IS is,PetscScalar *diag)
242: {
243: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
244: int i,ierr,N,*rows,m = A->m - 1,nz,*col;
245: PetscScalar *dd,*val;
248: ISGetLocalSize(is,&N);
249: ISGetIndices(is,&rows);
250: for (i=0; i<N; i++) {
251: if (rows[i]<0 || rows[i]>m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range");
252: MatGetRow(A,rows[i],&nz,&col,&val);
253: PetscMemzero(val,nz*sizeof(PetscScalar));
254: MatSetValues(A,1,&rows[i],nz,col,val,INSERT_VALUES);
255: MatRestoreRow(A,rows[i],&nz,&col,&val);
256: }
257: if (diag) {
258: if (a->mainbd == -1) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Main diagonal does not exist");
259: dd = a->diagv[a->mainbd];
260: for (i=0; i<N; i++) dd[rows[i]] = *diag;
261: }
262: ISRestoreIndices(is,&rows);
263: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
264: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
265: return(0);
266: }
268: #undef __FUNCT__
270: int MatGetSubMatrix_SeqBDiag(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *submat)
271: {
272: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
273: int nznew,*smap,i,j,ierr,oldcols = A->n;
274: int *irow,*icol,newr,newc,*cwork,*col,nz,bs;
275: PetscScalar *vwork,*val;
276: Mat newmat;
279: if (scall == MAT_REUSE_MATRIX) { /* no support for reuse so simply destroy all */
280: MatDestroy(*submat);
281: }
283: ISGetIndices(isrow,&irow);
284: ISGetIndices(iscol,&icol);
285: ISGetLocalSize(isrow,&newr);
286: ISGetLocalSize(iscol,&newc);
288: PetscMalloc((oldcols+1)*sizeof(int),&smap);
289: PetscMalloc((newc+1)*sizeof(int),&cwork);
290: PetscMalloc((newc+1)*sizeof(PetscScalar),&vwork);
291: ierr = PetscMemzero((char*)smap,oldcols*sizeof(int));
292: for (i=0; i<newc; i++) smap[icol[i]] = i+1;
294: /* Determine diagonals; then create submatrix */
295: bs = a->bs; /* Default block size remains the same */
296: MatCreateSeqBDiag(A->comm,newr,newc,0,bs,0,0,&newmat);
298: /* Fill new matrix */
299: for (i=0; i<newr; i++) {
300: MatGetRow(A,irow[i],&nz,&col,&val);
301: nznew = 0;
302: for (j=0; j<nz; j++) {
303: if (smap[col[j]]) {
304: cwork[nznew] = smap[col[j]] - 1;
305: vwork[nznew++] = val[j];
306: }
307: }
308: MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES);
309: MatRestoreRow(A,i,&nz,&col,&val);
310: }
311: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
312: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
314: /* Free work space */
315: PetscFree(smap);
316: PetscFree(cwork);
317: PetscFree(vwork);
318: ISRestoreIndices(isrow,&irow);
319: ISRestoreIndices(iscol,&icol);
320: *submat = newmat;
321: return(0);
322: }
324: #undef __FUNCT__
326: int MatGetSubMatrices_SeqBDiag(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B)
327: {
328: int ierr,i;
331: if (scall == MAT_INITIAL_MATRIX) {
332: PetscMalloc((n+1)*sizeof(Mat),B);
333: }
335: for (i=0; i<n; i++) {
336: MatGetSubMatrix_SeqBDiag(A,irow[i],icol[i],scall,&(*B)[i]);
337: }
338: return(0);
339: }
341: #undef __FUNCT__
343: int MatScale_SeqBDiag(PetscScalar *alpha,Mat inA)
344: {
345: Mat_SeqBDiag *a = (Mat_SeqBDiag*)inA->data;
346: int one = 1,i,len,bs = a->bs;
349: for (i=0; i<a->nd; i++) {
350: len = bs*bs*a->bdlen[i];
351: if (a->diag[i] > 0) {
352: BLscal_(&len,alpha,a->diagv[i] + bs*bs*a->diag[i],&one);
353: } else {
354: BLscal_(&len,alpha,a->diagv[i],&one);
355: }
356: }
357: PetscLogFlops(a->nz);
358: return(0);
359: }
361: #undef __FUNCT__
363: int MatDiagonalScale_SeqBDiag(Mat A,Vec ll,Vec rr)
364: {
365: Mat_SeqBDiag *a = (Mat_SeqBDiag*)A->data;
366: PetscScalar *l,*r,*dv;
367: int d,j,len,ierr;
368: int nd = a->nd,bs = a->bs,diag,m,n;
371: if (ll) {
372: VecGetSize(ll,&m);
373: if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length");
374: if (bs == 1) {
375: VecGetArray(ll,&l);
376: for (d=0; d<nd; d++) {
377: dv = a->diagv[d];
378: diag = a->diag[d];
379: len = a->bdlen[d];
380: if (diag > 0) for (j=0; j<len; j++) dv[j+diag] *= l[j+diag];
381: else for (j=0; j<len; j++) dv[j] *= l[j];
382: }
383: VecRestoreArray(ll,&l);
384: PetscLogFlops(a->nz);
385: } else SETERRQ(PETSC_ERR_SUP,"Not yet done for bs>1");
386: }
387: if (rr) {
388: VecGetSize(rr,&n);
389: if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length");
390: if (bs == 1) {
391: VecGetArray(rr,&r);
392: for (d=0; d<nd; d++) {
393: dv = a->diagv[d];
394: diag = a->diag[d];
395: len = a->bdlen[d];
396: if (diag > 0) for (j=0; j<len; j++) dv[j+diag] *= r[j];
397: else for (j=0; j<len; j++) dv[j] *= r[j-diag];
398: }
399: VecRestoreArray(rr,&r);
400: PetscLogFlops(a->nz);
401: } else SETERRQ(PETSC_ERR_SUP,"Not yet done for bs>1");
402: }
403: return(0);
404: }
406: static int MatDuplicate_SeqBDiag(Mat,MatDuplicateOption,Mat *);
407: EXTERN int MatLUFactorSymbolic_SeqBDiag(Mat,IS,IS,MatLUInfo*,Mat*);
408: EXTERN int MatILUFactorSymbolic_SeqBDiag(Mat,IS,IS,MatILUInfo*,Mat*);
409: EXTERN int MatILUFactor_SeqBDiag(Mat,IS,IS,MatILUInfo*);
410: EXTERN int MatLUFactorNumeric_SeqBDiag_N(Mat,Mat*);
411: EXTERN int MatLUFactorNumeric_SeqBDiag_1(Mat,Mat*);
412: EXTERN int MatSolve_SeqBDiag_1(Mat,Vec,Vec);
413: EXTERN int MatSolve_SeqBDiag_2(Mat,Vec,Vec);
414: EXTERN int MatSolve_SeqBDiag_3(Mat,Vec,Vec);
415: EXTERN int MatSolve_SeqBDiag_4(Mat,Vec,Vec);
416: EXTERN int MatSolve_SeqBDiag_5(Mat,Vec,Vec);
417: EXTERN int MatSolve_SeqBDiag_N(Mat,Vec,Vec);
419: #undef __FUNCT__
421: int MatSetUpPreallocation_SeqBDiag(Mat A)
422: {
423: int ierr;
426: MatSeqBDiagSetPreallocation(A,PETSC_DEFAULT,PETSC_DEFAULT,0,0);
427: return(0);
428: }
430: /* -------------------------------------------------------------------*/
431: static struct _MatOps MatOps_Values = {MatSetValues_SeqBDiag_N,
432: MatGetRow_SeqBDiag,
433: MatRestoreRow_SeqBDiag,
434: MatMult_SeqBDiag_N,
435: MatMultAdd_SeqBDiag_N,
436: MatMultTranspose_SeqBDiag_N,
437: MatMultTransposeAdd_SeqBDiag_N,
438: MatSolve_SeqBDiag_N,
439: 0,
440: 0,
441: 0,
442: 0,
443: 0,
444: MatRelax_SeqBDiag_N,
445: MatTranspose_SeqBDiag,
446: MatGetInfo_SeqBDiag,
447: 0,
448: MatGetDiagonal_SeqBDiag_N,
449: MatDiagonalScale_SeqBDiag,
450: MatNorm_SeqBDiag,
451: 0,
452: MatAssemblyEnd_SeqBDiag,
453: 0,
454: MatSetOption_SeqBDiag,
455: MatZeroEntries_SeqBDiag,
456: MatZeroRows_SeqBDiag,
457: 0,
458: MatLUFactorNumeric_SeqBDiag_N,
459: 0,
460: 0,
461: MatSetUpPreallocation_SeqBDiag,
462: MatILUFactorSymbolic_SeqBDiag,
463: 0,
464: 0,
465: 0,
466: MatDuplicate_SeqBDiag,
467: 0,
468: 0,
469: MatILUFactor_SeqBDiag,
470: 0,
471: 0,
472: MatGetSubMatrices_SeqBDiag,
473: 0,
474: MatGetValues_SeqBDiag_N,
475: 0,
476: MatPrintHelp_SeqBDiag,
477: MatScale_SeqBDiag,
478: 0,
479: 0,
480: 0,
481: MatGetBlockSize_SeqBDiag,
482: 0,
483: 0,
484: 0,
485: 0,
486: 0,
487: 0,
488: 0,
489: 0,
490: 0,
491: 0,
492: MatDestroy_SeqBDiag,
493: MatView_SeqBDiag,
494: MatGetPetscMaps_Petsc};
496: #undef __FUNCT__
498: /*@C
499: MatSeqBDiagSetPreallocation - Sets the nonzero structure and (optionally) arrays.
501: Collective on MPI_Comm
503: Input Parameters:
504: + B - the matrix
505: . nd - number of block diagonals (optional)
506: . bs - each element of a diagonal is an bs x bs dense matrix
507: . diag - optional array of block diagonal numbers (length nd).
508: For a matrix element A[i,j], where i=row and j=column, the
509: diagonal number is
510: $ diag = i/bs - j/bs (integer division)
511: Set diag=PETSC_NULL on input for PETSc to dynamically allocate memory as
512: needed (expensive).
513: - diagv - pointer to actual diagonals (in same order as diag array),
514: if allocated by user. Otherwise, set diagv=PETSC_NULL on input for PETSc
515: to control memory allocation.
517: Options Database Keys:
518: . -mat_block_size <bs> - Sets blocksize
519: . -mat_bdiag_diags <s1,s2,s3,...> - Sets diagonal numbers
521: Notes:
522: See the users manual for further details regarding this storage format.
524: Fortran Note:
525: Fortran programmers cannot set diagv; this value is ignored.
527: Level: intermediate
529: .keywords: matrix, block, diagonal, sparse
531: .seealso: MatCreate(), MatCreateMPIBDiag(), MatSetValues()
532: @*/
533: int MatSeqBDiagSetPreallocation(Mat B,int nd,int bs,int *diag,PetscScalar **diagv)
534: {
535: Mat_SeqBDiag *b;
536: int i,nda,sizetot,ierr, nd2 = 128,idiag[128];
537: PetscTruth flg1;
540: PetscTypeCompare((PetscObject)B,MATSEQBDIAG,&flg1);
541: if (!flg1) return(0);
543: B->preallocated = PETSC_TRUE;
544: if (bs == PETSC_DEFAULT) bs = 1;
545: if (bs == 0) SETERRQ(1,"Blocksize cannot be zero");
546: if (nd == PETSC_DEFAULT) nd = 0;
547: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
548: PetscOptionsGetIntArray(PETSC_NULL,"-mat_bdiag_diags",idiag,&nd2,&flg1);
549: if (flg1) {
550: diag = idiag;
551: nd = nd2;
552: }
554: if ((B->n%bs) || (B->m%bs)) SETERRQ(PETSC_ERR_ARG_SIZ,"Invalid block size");
555: if (!nd) nda = nd + 1;
556: else nda = nd;
557: b = (Mat_SeqBDiag*)B->data;
559: PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg1);
560: if (!flg1) {
561: switch (bs) {
562: case 1:
563: B->ops->setvalues = MatSetValues_SeqBDiag_1;
564: B->ops->getvalues = MatGetValues_SeqBDiag_1;
565: B->ops->getdiagonal = MatGetDiagonal_SeqBDiag_1;
566: B->ops->mult = MatMult_SeqBDiag_1;
567: B->ops->multadd = MatMultAdd_SeqBDiag_1;
568: B->ops->multtranspose = MatMultTranspose_SeqBDiag_1;
569: B->ops->multtransposeadd= MatMultTransposeAdd_SeqBDiag_1;
570: B->ops->relax = MatRelax_SeqBDiag_1;
571: B->ops->solve = MatSolve_SeqBDiag_1;
572: B->ops->lufactornumeric = MatLUFactorNumeric_SeqBDiag_1;
573: break;
574: case 2:
575: B->ops->mult = MatMult_SeqBDiag_2;
576: B->ops->multadd = MatMultAdd_SeqBDiag_2;
577: B->ops->solve = MatSolve_SeqBDiag_2;
578: break;
579: case 3:
580: B->ops->mult = MatMult_SeqBDiag_3;
581: B->ops->multadd = MatMultAdd_SeqBDiag_3;
582: B->ops->solve = MatSolve_SeqBDiag_3;
583: break;
584: case 4:
585: B->ops->mult = MatMult_SeqBDiag_4;
586: B->ops->multadd = MatMultAdd_SeqBDiag_4;
587: B->ops->solve = MatSolve_SeqBDiag_4;
588: break;
589: case 5:
590: B->ops->mult = MatMult_SeqBDiag_5;
591: B->ops->multadd = MatMultAdd_SeqBDiag_5;
592: B->ops->solve = MatSolve_SeqBDiag_5;
593: break;
594: }
595: }
597: b->mblock = B->m/bs;
598: b->nblock = B->n/bs;
599: b->nd = nd;
600: b->bs = bs;
601: b->ndim = 0;
602: b->mainbd = -1;
603: b->pivot = 0;
605: ierr = PetscMalloc(2*nda*sizeof(int),&b->diag);
606: b->bdlen = b->diag + nda;
607: ierr = PetscMalloc((B->n+1)*sizeof(int),&b->colloc);
608: ierr = PetscMalloc(nda*sizeof(PetscScalar*),&b->diagv);
609: sizetot = 0;
611: if (diagv) { /* user allocated space */
612: b->user_alloc = PETSC_TRUE;
613: for (i=0; i<nd; i++) b->diagv[i] = diagv[i];
614: } else b->user_alloc = PETSC_FALSE;
616: for (i=0; i<nd; i++) {
617: b->diag[i] = diag[i];
618: if (diag[i] > 0) { /* lower triangular */
619: b->bdlen[i] = PetscMin(b->nblock,b->mblock - diag[i]);
620: } else { /* upper triangular */
621: b->bdlen[i] = PetscMin(b->mblock,b->nblock + diag[i]);
622: }
623: sizetot += b->bdlen[i];
624: }
625: sizetot *= bs*bs;
626: b->maxnz = sizetot;
627: ierr = PetscMalloc((B->n+1)*sizeof(PetscScalar),&b->dvalue);
628: PetscLogObjectMemory(B,(nda*(bs+2))*sizeof(int) + bs*nda*sizeof(PetscScalar)
629: + nda*sizeof(PetscScalar*) + sizeof(Mat_SeqBDiag)
630: + sizeof(struct _p_Mat) + sizetot*sizeof(PetscScalar));
632: if (!b->user_alloc) {
633: for (i=0; i<nd; i++) {
634: PetscMalloc(bs*bs*b->bdlen[i]*sizeof(PetscScalar),&b->diagv[i]);
635: PetscMemzero(b->diagv[i],bs*bs*b->bdlen[i]*sizeof(PetscScalar));
636: }
637: b->nonew = 0; b->nonew_diag = 0;
638: } else { /* diagonals are set on input; don't allow dynamic allocation */
639: b->nonew = 1; b->nonew_diag = 1;
640: }
642: /* adjust diagv so one may access rows with diagv[diag][row] for all rows */
643: for (i=0; i<nd; i++) {
644: if (diag[i] > 0) {
645: b->diagv[i] -= bs*bs*diag[i];
646: }
647: }
649: b->nz = b->maxnz; /* Currently not keeping track of exact count */
650: b->roworiented = PETSC_TRUE;
651: B->info.nz_unneeded = (double)b->maxnz;
652: return(0);
653: }
655: #undef __FUNCT__
657: static int MatDuplicate_SeqBDiag(Mat A,MatDuplicateOption cpvalues,Mat *matout)
658: {
659: Mat_SeqBDiag *newmat,*a = (Mat_SeqBDiag*)A->data;
660: int i,ierr,len,diag,bs = a->bs;
661: Mat mat;
664: MatCreateSeqBDiag(A->comm,A->m,A->n,a->nd,bs,a->diag,PETSC_NULL,matout);
666: /* Copy contents of diagonals */
667: mat = *matout;
668: newmat = (Mat_SeqBDiag*)mat->data;
669: if (cpvalues == MAT_COPY_VALUES) {
670: for (i=0; i<a->nd; i++) {
671: len = a->bdlen[i] * bs * bs * sizeof(PetscScalar);
672: diag = a->diag[i];
673: if (diag > 0) {
674: PetscMemcpy(newmat->diagv[i]+bs*bs*diag,a->diagv[i]+bs*bs*diag,len);
675: } else {
676: PetscMemcpy(newmat->diagv[i],a->diagv[i],len);
677: }
678: }
679: }
680: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
681: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
682: return(0);
683: }
685: EXTERN_C_BEGIN
686: #undef __FUNCT__
688: int MatLoad_SeqBDiag(PetscViewer viewer,MatType type,Mat *A)
689: {
690: Mat B;
691: int *scols,i,nz,ierr,fd,header[4],size,nd = 128;
692: int bs,*rowlengths = 0,M,N,*cols,extra_rows,*diag = 0;
693: int idiag[128];
694: PetscScalar *vals,*svals;
695: MPI_Comm comm;
696: PetscTruth flg;
697:
699: PetscObjectGetComm((PetscObject)viewer,&comm);
700: MPI_Comm_size(comm,&size);
701: if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor");
702: PetscViewerBinaryGetDescriptor(viewer,&fd);
703: PetscBinaryRead(fd,header,4,PETSC_INT);
704: if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
705: M = header[1]; N = header[2]; nz = header[3];
706: if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only load square matrices");
707: if (header[3] < 0) {
708: SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBDiag");
709: }
711: /*
712: This code adds extra rows to make sure the number of rows is
713: divisible by the blocksize
714: */
715: bs = 1;
716: PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);
717: extra_rows = bs - M + bs*(M/bs);
718: if (extra_rows == bs) extra_rows = 0;
719: if (extra_rows) {
720: PetscLogInfo(0,"MatLoad_SeqBDiag:Padding loaded matrix to match blocksizen");
721: }
723: /* read row lengths */
724: PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);
725: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
726: for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;
728: /* load information about diagonals */
729: PetscOptionsGetIntArray(PETSC_NULL,"-matload_bdiag_diags",idiag,&nd,&flg);
730: if (flg) {
731: diag = idiag;
732: }
734: /* create our matrix */
735: MatCreateSeqBDiag(comm,M+extra_rows,M+extra_rows,nd,bs,diag,
736: PETSC_NULL,A);
737: B = *A;
739: /* read column indices and nonzeros */
740: PetscMalloc(nz*sizeof(int),&scols);
741: cols = scols;
742: PetscBinaryRead(fd,cols,nz,PETSC_INT);
743: PetscMalloc(nz*sizeof(PetscScalar),&svals);
744: vals = svals;
745: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
746: /* insert into matrix */
748: for (i=0; i<M; i++) {
749: MatSetValues(B,1,&i,rowlengths[i],scols,svals,INSERT_VALUES);
750: scols += rowlengths[i]; svals += rowlengths[i];
751: }
752: vals[0] = 1.0;
753: for (i=M; i<M+extra_rows; i++) {
754: MatSetValues(B,1,&i,1,&i,vals,INSERT_VALUES);
755: }
757: PetscFree(cols);
758: PetscFree(vals);
759: PetscFree(rowlengths);
761: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
762: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
763: return(0);
764: }
765: EXTERN_C_END
767: EXTERN_C_BEGIN
768: #undef __FUNCT__
770: int MatCreate_SeqBDiag(Mat B)
771: {
772: Mat_SeqBDiag *b;
773: int ierr,size;
776: MPI_Comm_size(B->comm,&size);
777: if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
779: B->m = B->M = PetscMax(B->m,B->M);
780: B->n = B->N = PetscMax(B->n,B->N);
782: ierr = PetscNew(Mat_SeqBDiag,&b);
783: B->data = (void*)b;
784: ierr = PetscMemzero(b,sizeof(Mat_SeqBDiag));
785: ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
786: B->factor = 0;
787: B->mapping = 0;
789: PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);
790: PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);
792: b->ndim = 0;
793: b->mainbd = -1;
794: b->pivot = 0;
796: b->roworiented = PETSC_TRUE;
797: return(0);
798: }
799: EXTERN_C_END
801: #undef __FUNCT__
803: /*@C
804: MatCreateSeqBDiag - Creates a sequential block diagonal matrix.
806: Collective on MPI_Comm
808: Input Parameters:
809: + comm - MPI communicator, set to PETSC_COMM_SELF
810: . m - number of rows
811: . n - number of columns
812: . nd - number of block diagonals (optional)
813: . bs - each element of a diagonal is an bs x bs dense matrix
814: . diag - optional array of block diagonal numbers (length nd).
815: For a matrix element A[i,j], where i=row and j=column, the
816: diagonal number is
817: $ diag = i/bs - j/bs (integer division)
818: Set diag=PETSC_NULL on input for PETSc to dynamically allocate memory as
819: needed (expensive).
820: - diagv - pointer to actual diagonals (in same order as diag array),
821: if allocated by user. Otherwise, set diagv=PETSC_NULL on input for PETSc
822: to control memory allocation.
824: Output Parameters:
825: . A - the matrix
827: Options Database Keys:
828: . -mat_block_size <bs> - Sets blocksize
829: . -mat_bdiag_diags <s1,s2,s3,...> - Sets diagonal numbers
831: Notes:
832: See the users manual for further details regarding this storage format.
834: Fortran Note:
835: Fortran programmers cannot set diagv; this value is ignored.
837: Level: intermediate
839: .keywords: matrix, block, diagonal, sparse
841: .seealso: MatCreate(), MatCreateMPIBDiag(), MatSetValues()
842: @*/
843: int MatCreateSeqBDiag(MPI_Comm comm,int m,int n,int nd,int bs,int *diag,PetscScalar **diagv,Mat *A)
844: {
848: MatCreate(comm,m,n,m,n,A);
849: MatSetType(*A,MATSEQBDIAG);
850: MatSeqBDiagSetPreallocation(*A,nd,bs,diag,diagv);
851: return(0);
852: }