Actual source code: aij.c
1: /*
2: Defines the basic matrix operations for the AIJ (compressed row)
3: matrix storage format.
4: */
6: #include <../src/mat/impls/aij/seq/aij.h>
7: #include <petscblaslapack.h>
8: #include <petscbt.h>
9: #include <petsc/private/kernels/blocktranspose.h>
11: PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A)
12: {
13: PetscErrorCode ierr;
14: PetscBool flg;
15: char type[256];
17: PetscObjectOptionsBegin((PetscObject)A);
18: PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg);
19: if (flg) MatSeqAIJSetType(A,type);
20: PetscOptionsEnd();
21: return 0;
22: }
24: PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions)
25: {
26: PetscInt i,m,n;
27: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
29: MatGetSize(A,&m,&n);
30: PetscArrayzero(reductions,n);
31: if (type == NORM_2) {
32: for (i=0; i<aij->i[m]; i++) {
33: reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]);
34: }
35: } else if (type == NORM_1) {
36: for (i=0; i<aij->i[m]; i++) {
37: reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]);
38: }
39: } else if (type == NORM_INFINITY) {
40: for (i=0; i<aij->i[m]; i++) {
41: reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]);
42: }
43: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
44: for (i=0; i<aij->i[m]; i++) {
45: reductions[aij->j[i]] += PetscRealPart(aij->a[i]);
46: }
47: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
48: for (i=0; i<aij->i[m]; i++) {
49: reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]);
50: }
51: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type");
53: if (type == NORM_2) {
54: for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
55: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
56: for (i=0; i<n; i++) reductions[i] /= m;
57: }
58: return 0;
59: }
61: PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is)
62: {
63: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
64: PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs;
65: const PetscInt *jj = a->j,*ii = a->i;
66: PetscInt *rows;
68: for (i=0; i<m; i++) {
69: if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) {
70: cnt++;
71: }
72: }
73: PetscMalloc1(cnt,&rows);
74: cnt = 0;
75: for (i=0; i<m; i++) {
76: if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) {
77: rows[cnt] = i;
78: cnt++;
79: }
80: }
81: ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);
82: return 0;
83: }
85: PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows)
86: {
87: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
88: const MatScalar *aa;
89: PetscInt i,m=A->rmap->n,cnt = 0;
90: const PetscInt *ii = a->i,*jj = a->j,*diag;
91: PetscInt *rows;
93: MatSeqAIJGetArrayRead(A,&aa);
94: MatMarkDiagonal_SeqAIJ(A);
95: diag = a->diag;
96: for (i=0; i<m; i++) {
97: if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
98: cnt++;
99: }
100: }
101: PetscMalloc1(cnt,&rows);
102: cnt = 0;
103: for (i=0; i<m; i++) {
104: if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
105: rows[cnt++] = i;
106: }
107: }
108: *nrows = cnt;
109: *zrows = rows;
110: MatSeqAIJRestoreArrayRead(A,&aa);
111: return 0;
112: }
114: PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows)
115: {
116: PetscInt nrows,*rows;
118: *zrows = NULL;
119: MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);
120: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);
121: return 0;
122: }
124: PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows)
125: {
126: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
127: const MatScalar *aa;
128: PetscInt m=A->rmap->n,cnt = 0;
129: const PetscInt *ii;
130: PetscInt n,i,j,*rows;
132: MatSeqAIJGetArrayRead(A,&aa);
133: *keptrows = NULL;
134: ii = a->i;
135: for (i=0; i<m; i++) {
136: n = ii[i+1] - ii[i];
137: if (!n) {
138: cnt++;
139: goto ok1;
140: }
141: for (j=ii[i]; j<ii[i+1]; j++) {
142: if (aa[j] != 0.0) goto ok1;
143: }
144: cnt++;
145: ok1:;
146: }
147: if (!cnt) {
148: MatSeqAIJRestoreArrayRead(A,&aa);
149: return 0;
150: }
151: PetscMalloc1(A->rmap->n-cnt,&rows);
152: cnt = 0;
153: for (i=0; i<m; i++) {
154: n = ii[i+1] - ii[i];
155: if (!n) continue;
156: for (j=ii[i]; j<ii[i+1]; j++) {
157: if (aa[j] != 0.0) {
158: rows[cnt++] = i;
159: break;
160: }
161: }
162: }
163: MatSeqAIJRestoreArrayRead(A,&aa);
164: ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);
165: return 0;
166: }
168: PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is)
169: {
170: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data;
171: PetscInt i,m = Y->rmap->n;
172: const PetscInt *diag;
173: MatScalar *aa;
174: const PetscScalar *v;
175: PetscBool missing;
177: if (Y->assembled) {
178: MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);
179: if (!missing) {
180: diag = aij->diag;
181: VecGetArrayRead(D,&v);
182: MatSeqAIJGetArray(Y,&aa);
183: if (is == INSERT_VALUES) {
184: for (i=0; i<m; i++) {
185: aa[diag[i]] = v[i];
186: }
187: } else {
188: for (i=0; i<m; i++) {
189: aa[diag[i]] += v[i];
190: }
191: }
192: MatSeqAIJRestoreArray(Y,&aa);
193: VecRestoreArrayRead(D,&v);
194: return 0;
195: }
196: MatSeqAIJInvalidateDiagonal(Y);
197: }
198: MatDiagonalSet_Default(Y,D,is);
199: return 0;
200: }
202: PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
203: {
204: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
205: PetscInt i,ishift;
207: *m = A->rmap->n;
208: if (!ia) return 0;
209: ishift = 0;
210: if (symmetric && !A->structurally_symmetric) {
211: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);
212: } else if (oshift == 1) {
213: PetscInt *tia;
214: PetscInt nz = a->i[A->rmap->n];
215: /* malloc space and add 1 to i and j indices */
216: PetscMalloc1(A->rmap->n+1,&tia);
217: for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1;
218: *ia = tia;
219: if (ja) {
220: PetscInt *tja;
221: PetscMalloc1(nz+1,&tja);
222: for (i=0; i<nz; i++) tja[i] = a->j[i] + 1;
223: *ja = tja;
224: }
225: } else {
226: *ia = a->i;
227: if (ja) *ja = a->j;
228: }
229: return 0;
230: }
232: PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
233: {
234: if (!ia) return 0;
235: if ((symmetric && !A->structurally_symmetric) || oshift == 1) {
236: PetscFree(*ia);
237: if (ja) PetscFree(*ja);
238: }
239: return 0;
240: }
242: PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
243: {
244: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
245: PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n;
246: PetscInt nz = a->i[m],row,*jj,mr,col;
248: *nn = n;
249: if (!ia) return 0;
250: if (symmetric) {
251: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja);
252: } else {
253: PetscCalloc1(n,&collengths);
254: PetscMalloc1(n+1,&cia);
255: PetscMalloc1(nz,&cja);
256: jj = a->j;
257: for (i=0; i<nz; i++) {
258: collengths[jj[i]]++;
259: }
260: cia[0] = oshift;
261: for (i=0; i<n; i++) {
262: cia[i+1] = cia[i] + collengths[i];
263: }
264: PetscArrayzero(collengths,n);
265: jj = a->j;
266: for (row=0; row<m; row++) {
267: mr = a->i[row+1] - a->i[row];
268: for (i=0; i<mr; i++) {
269: col = *jj++;
271: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
272: }
273: }
274: PetscFree(collengths);
275: *ia = cia; *ja = cja;
276: }
277: return 0;
278: }
280: PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
281: {
282: if (!ia) return 0;
284: PetscFree(*ia);
285: PetscFree(*ja);
286: return 0;
287: }
289: /*
290: MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from
291: MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output
292: spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ()
293: */
294: PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
295: {
296: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
297: PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n;
298: PetscInt nz = a->i[m],row,mr,col,tmp;
299: PetscInt *cspidx;
300: const PetscInt *jj;
302: *nn = n;
303: if (!ia) return 0;
305: PetscCalloc1(n,&collengths);
306: PetscMalloc1(n+1,&cia);
307: PetscMalloc1(nz,&cja);
308: PetscMalloc1(nz,&cspidx);
309: jj = a->j;
310: for (i=0; i<nz; i++) {
311: collengths[jj[i]]++;
312: }
313: cia[0] = oshift;
314: for (i=0; i<n; i++) {
315: cia[i+1] = cia[i] + collengths[i];
316: }
317: PetscArrayzero(collengths,n);
318: jj = a->j;
319: for (row=0; row<m; row++) {
320: mr = a->i[row+1] - a->i[row];
321: for (i=0; i<mr; i++) {
322: col = *jj++;
323: tmp = cia[col] + collengths[col]++ - oshift;
324: cspidx[tmp] = a->i[row] + i; /* index of a->j */
325: cja[tmp] = row + oshift;
326: }
327: }
328: PetscFree(collengths);
329: *ia = cia;
330: *ja = cja;
331: *spidx = cspidx;
332: return 0;
333: }
335: PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
336: {
337: MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);
338: PetscFree(*spidx);
339: return 0;
340: }
342: PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[])
343: {
344: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
345: PetscInt *ai = a->i;
346: PetscScalar *aa;
348: MatSeqAIJGetArray(A,&aa);
349: PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row]);
350: MatSeqAIJRestoreArray(A,&aa);
351: return 0;
352: }
354: /*
355: MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions
357: - a single row of values is set with each call
358: - no row or column indices are negative or (in error) larger than the number of rows or columns
359: - the values are always added to the matrix, not set
360: - no new locations are introduced in the nonzero structure of the matrix
362: This does NOT assume the global column indices are sorted
364: */
366: #include <petsc/private/isimpl.h>
367: PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
368: {
369: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
370: PetscInt low,high,t,row,nrow,i,col,l;
371: const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j;
372: PetscInt lastcol = -1;
373: MatScalar *ap,value,*aa;
374: const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices;
376: MatSeqAIJGetArray(A,&aa);
377: row = ridx[im[0]];
378: rp = aj + ai[row];
379: ap = aa + ai[row];
380: nrow = ailen[row];
381: low = 0;
382: high = nrow;
383: for (l=0; l<n; l++) { /* loop over added columns */
384: col = cidx[in[l]];
385: value = v[l];
387: if (col <= lastcol) low = 0;
388: else high = nrow;
389: lastcol = col;
390: while (high-low > 5) {
391: t = (low+high)/2;
392: if (rp[t] > col) high = t;
393: else low = t;
394: }
395: for (i=low; i<high; i++) {
396: if (rp[i] == col) {
397: ap[i] += value;
398: low = i + 1;
399: break;
400: }
401: }
402: }
403: MatSeqAIJRestoreArray(A,&aa);
404: return 0;
405: }
407: PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
408: {
409: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
410: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
411: PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen;
412: PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1;
413: MatScalar *ap=NULL,value=0.0,*aa;
414: PetscBool ignorezeroentries = a->ignorezeroentries;
415: PetscBool roworiented = a->roworiented;
417: MatSeqAIJGetArray(A,&aa);
418: for (k=0; k<m; k++) { /* loop over added rows */
419: row = im[k];
420: if (row < 0) continue;
422: rp = aj + ai[row];
423: if (!A->structure_only) ap = aa + ai[row];
424: rmax = imax[row]; nrow = ailen[row];
425: low = 0;
426: high = nrow;
427: for (l=0; l<n; l++) { /* loop over added columns */
428: if (in[l] < 0) continue;
430: col = in[l];
431: if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m];
432: if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue;
434: if (col <= lastcol) low = 0;
435: else high = nrow;
436: lastcol = col;
437: while (high-low > 5) {
438: t = (low+high)/2;
439: if (rp[t] > col) high = t;
440: else low = t;
441: }
442: for (i=low; i<high; i++) {
443: if (rp[i] > col) break;
444: if (rp[i] == col) {
445: if (!A->structure_only) {
446: if (is == ADD_VALUES) {
447: ap[i] += value;
448: (void)PetscLogFlops(1.0);
449: }
450: else ap[i] = value;
451: }
452: low = i + 1;
453: goto noinsert;
454: }
455: }
456: if (value == 0.0 && ignorezeroentries && row != col) goto noinsert;
457: if (nonew == 1) goto noinsert;
459: if (A->structure_only) {
460: MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar);
461: } else {
462: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
463: }
464: N = nrow++ - 1; a->nz++; high++;
465: /* shift up all the later entries in this row */
466: PetscArraymove(rp+i+1,rp+i,N-i+1);
467: rp[i] = col;
468: if (!A->structure_only) {
469: PetscArraymove(ap+i+1,ap+i,N-i+1);
470: ap[i] = value;
471: }
472: low = i + 1;
473: A->nonzerostate++;
474: noinsert:;
475: }
476: ailen[row] = nrow;
477: }
478: MatSeqAIJRestoreArray(A,&aa);
479: return 0;
480: }
482: PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
483: {
484: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
485: PetscInt *rp,k,row;
486: PetscInt *ai = a->i;
487: PetscInt *aj = a->j;
488: MatScalar *aa,*ap;
493: MatSeqAIJGetArray(A,&aa);
494: for (k=0; k<m; k++) { /* loop over added rows */
495: row = im[k];
496: rp = aj + ai[row];
497: ap = aa + ai[row];
499: PetscMemcpy(rp,in,n*sizeof(PetscInt));
500: if (!A->structure_only) {
501: if (v) {
502: PetscMemcpy(ap,v,n*sizeof(PetscScalar));
503: v += n;
504: } else {
505: PetscMemzero(ap,n*sizeof(PetscScalar));
506: }
507: }
508: a->ilen[row] = n;
509: a->imax[row] = n;
510: a->i[row+1] = a->i[row]+n;
511: a->nz += n;
512: }
513: MatSeqAIJRestoreArray(A,&aa);
514: return 0;
515: }
517: /*@
518: MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix.
520: Input Parameters:
521: + A - the SeqAIJ matrix
522: - nztotal - bound on the number of nonzeros
524: Level: advanced
526: Notes:
527: This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row.
528: Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used
529: as always with multiple matrix assemblies.
531: .seealso: MatSetOption(), MAT_SORTED_FULL, MatSetValues(), MatSeqAIJSetPreallocation()
532: @*/
534: PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal)
535: {
536: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
538: PetscLayoutSetUp(A->rmap);
539: PetscLayoutSetUp(A->cmap);
540: a->maxnz = nztotal;
541: if (!a->imax) {
542: PetscMalloc1(A->rmap->n,&a->imax);
543: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt));
544: }
545: if (!a->ilen) {
546: PetscMalloc1(A->rmap->n,&a->ilen);
547: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt));
548: } else {
549: PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));
550: }
552: /* allocate the matrix space */
553: if (A->structure_only) {
554: PetscMalloc1(nztotal,&a->j);
555: PetscMalloc1(A->rmap->n+1,&a->i);
556: PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt));
557: } else {
558: PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i);
559: PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)));
560: }
561: a->i[0] = 0;
562: if (A->structure_only) {
563: a->singlemalloc = PETSC_FALSE;
564: a->free_a = PETSC_FALSE;
565: } else {
566: a->singlemalloc = PETSC_TRUE;
567: a->free_a = PETSC_TRUE;
568: }
569: a->free_ij = PETSC_TRUE;
570: A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation;
571: A->preallocated = PETSC_TRUE;
572: return 0;
573: }
575: PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
576: {
577: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
578: PetscInt *rp,k,row;
579: PetscInt *ai = a->i,*ailen = a->ilen;
580: PetscInt *aj = a->j;
581: MatScalar *aa,*ap;
583: MatSeqAIJGetArray(A,&aa);
584: for (k=0; k<m; k++) { /* loop over added rows */
585: row = im[k];
587: rp = aj + ai[row];
588: ap = aa + ai[row];
589: if (!A->was_assembled) {
590: PetscMemcpy(rp,in,n*sizeof(PetscInt));
591: }
592: if (!A->structure_only) {
593: if (v) {
594: PetscMemcpy(ap,v,n*sizeof(PetscScalar));
595: v += n;
596: } else {
597: PetscMemzero(ap,n*sizeof(PetscScalar));
598: }
599: }
600: ailen[row] = n;
601: a->nz += n;
602: }
603: MatSeqAIJRestoreArray(A,&aa);
604: return 0;
605: }
607: PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
608: {
609: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
610: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
611: PetscInt *ai = a->i,*ailen = a->ilen;
612: MatScalar *ap,*aa;
614: MatSeqAIJGetArray(A,&aa);
615: for (k=0; k<m; k++) { /* loop over rows */
616: row = im[k];
617: if (row < 0) {v += n; continue;} /* negative row */
619: rp = aj + ai[row]; ap = aa + ai[row];
620: nrow = ailen[row];
621: for (l=0; l<n; l++) { /* loop over columns */
622: if (in[l] < 0) {v++; continue;} /* negative column */
624: col = in[l];
625: high = nrow; low = 0; /* assume unsorted */
626: while (high-low > 5) {
627: t = (low+high)/2;
628: if (rp[t] > col) high = t;
629: else low = t;
630: }
631: for (i=low; i<high; i++) {
632: if (rp[i] > col) break;
633: if (rp[i] == col) {
634: *v++ = ap[i];
635: goto finished;
636: }
637: }
638: *v++ = 0.0;
639: finished:;
640: }
641: }
642: MatSeqAIJRestoreArray(A,&aa);
643: return 0;
644: }
646: PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer)
647: {
648: Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data;
649: const PetscScalar *av;
650: PetscInt header[4],M,N,m,nz,i;
651: PetscInt *rowlens;
653: PetscViewerSetUp(viewer);
655: M = mat->rmap->N;
656: N = mat->cmap->N;
657: m = mat->rmap->n;
658: nz = A->nz;
660: /* write matrix header */
661: header[0] = MAT_FILE_CLASSID;
662: header[1] = M; header[2] = N; header[3] = nz;
663: PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);
665: /* fill in and store row lengths */
666: PetscMalloc1(m,&rowlens);
667: for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i];
668: PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT);
669: PetscFree(rowlens);
670: /* store column indices */
671: PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT);
672: /* store nonzero values */
673: MatSeqAIJGetArrayRead(mat,&av);
674: PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR);
675: MatSeqAIJRestoreArrayRead(mat,&av);
677: /* write block size option to the viewer's .info file */
678: MatView_Binary_BlockSizes(mat,viewer);
679: return 0;
680: }
682: static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer)
683: {
684: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
685: PetscInt i,k,m=A->rmap->N;
687: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
688: for (i=0; i<m; i++) {
689: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
690: for (k=a->i[i]; k<a->i[i+1]; k++) {
691: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k]);
692: }
693: PetscViewerASCIIPrintf(viewer,"\n");
694: }
695: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
696: return 0;
697: }
699: extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer);
701: PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer)
702: {
703: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
704: const PetscScalar *av;
705: PetscInt i,j,m = A->rmap->n;
706: const char *name;
707: PetscViewerFormat format;
709: if (A->structure_only) {
710: MatView_SeqAIJ_ASCII_structonly(A,viewer);
711: return 0;
712: }
714: PetscViewerGetFormat(viewer,&format);
715: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
717: /* trigger copy to CPU if needed */
718: MatSeqAIJGetArrayRead(A,&av);
719: MatSeqAIJRestoreArrayRead(A,&av);
720: if (format == PETSC_VIEWER_ASCII_MATLAB) {
721: PetscInt nofinalvalue = 0;
722: if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) {
723: /* Need a dummy value to ensure the dimension of the matrix. */
724: nofinalvalue = 1;
725: }
726: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
727: PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n);
728: PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz);
729: #if defined(PETSC_USE_COMPLEX)
730: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue);
731: #else
732: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue);
733: #endif
734: PetscViewerASCIIPrintf(viewer,"zzz = [\n");
736: for (i=0; i<m; i++) {
737: for (j=a->i[i]; j<a->i[i+1]; j++) {
738: #if defined(PETSC_USE_COMPLEX)
739: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
740: #else
741: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);
742: #endif
743: }
744: }
745: if (nofinalvalue) {
746: #if defined(PETSC_USE_COMPLEX)
747: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.);
748: #else
749: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0);
750: #endif
751: }
752: PetscObjectGetName((PetscObject)A,&name);
753: PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);
754: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
755: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
756: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
757: for (i=0; i<m; i++) {
758: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
759: for (j=a->i[i]; j<a->i[i+1]; j++) {
760: #if defined(PETSC_USE_COMPLEX)
761: if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) {
762: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
763: } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) {
764: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));
765: } else if (PetscRealPart(a->a[j]) != 0.0) {
766: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
767: }
768: #else
769: if (a->a[j] != 0.0) PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j]);
770: #endif
771: }
772: PetscViewerASCIIPrintf(viewer,"\n");
773: }
774: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
775: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
776: PetscInt nzd=0,fshift=1,*sptr;
777: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
778: PetscMalloc1(m+1,&sptr);
779: for (i=0; i<m; i++) {
780: sptr[i] = nzd+1;
781: for (j=a->i[i]; j<a->i[i+1]; j++) {
782: if (a->j[j] >= i) {
783: #if defined(PETSC_USE_COMPLEX)
784: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++;
785: #else
786: if (a->a[j] != 0.0) nzd++;
787: #endif
788: }
789: }
790: }
791: sptr[m] = nzd+1;
792: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd);
793: for (i=0; i<m+1; i+=6) {
794: if (i+4<m) {
795: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);
796: } else if (i+3<m) {
797: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);
798: } else if (i+2<m) {
799: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);
800: } else if (i+1<m) {
801: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2]);
802: } else if (i<m) {
803: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1]);
804: } else {
805: PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i]);
806: }
807: }
808: PetscViewerASCIIPrintf(viewer,"\n");
809: PetscFree(sptr);
810: for (i=0; i<m; i++) {
811: for (j=a->i[i]; j<a->i[i+1]; j++) {
812: if (a->j[j] >= i) PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift);
813: }
814: PetscViewerASCIIPrintf(viewer,"\n");
815: }
816: PetscViewerASCIIPrintf(viewer,"\n");
817: for (i=0; i<m; i++) {
818: for (j=a->i[i]; j<a->i[i+1]; j++) {
819: if (a->j[j] >= i) {
820: #if defined(PETSC_USE_COMPLEX)
821: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) {
822: PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
823: }
824: #else
825: if (a->a[j] != 0.0) PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);
826: #endif
827: }
828: }
829: PetscViewerASCIIPrintf(viewer,"\n");
830: }
831: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
832: } else if (format == PETSC_VIEWER_ASCII_DENSE) {
833: PetscInt cnt = 0,jcnt;
834: PetscScalar value;
835: #if defined(PETSC_USE_COMPLEX)
836: PetscBool realonly = PETSC_TRUE;
838: for (i=0; i<a->i[m]; i++) {
839: if (PetscImaginaryPart(a->a[i]) != 0.0) {
840: realonly = PETSC_FALSE;
841: break;
842: }
843: }
844: #endif
846: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
847: for (i=0; i<m; i++) {
848: jcnt = 0;
849: for (j=0; j<A->cmap->n; j++) {
850: if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) {
851: value = a->a[cnt++];
852: jcnt++;
853: } else {
854: value = 0.0;
855: }
856: #if defined(PETSC_USE_COMPLEX)
857: if (realonly) {
858: PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));
859: } else {
860: PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));
861: }
862: #else
863: PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);
864: #endif
865: }
866: PetscViewerASCIIPrintf(viewer,"\n");
867: }
868: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
869: } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) {
870: PetscInt fshift=1;
871: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
872: #if defined(PETSC_USE_COMPLEX)
873: PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n");
874: #else
875: PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n");
876: #endif
877: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz);
878: for (i=0; i<m; i++) {
879: for (j=a->i[i]; j<a->i[i+1]; j++) {
880: #if defined(PETSC_USE_COMPLEX)
881: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
882: #else
883: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);
884: #endif
885: }
886: }
887: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
888: } else {
889: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
890: if (A->factortype) {
891: for (i=0; i<m; i++) {
892: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
893: /* L part */
894: for (j=a->i[i]; j<a->i[i+1]; j++) {
895: #if defined(PETSC_USE_COMPLEX)
896: if (PetscImaginaryPart(a->a[j]) > 0.0) {
897: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
898: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
899: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));
900: } else {
901: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
902: }
903: #else
904: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j]);
905: #endif
906: }
907: /* diagonal */
908: j = a->diag[i];
909: #if defined(PETSC_USE_COMPLEX)
910: if (PetscImaginaryPart(a->a[j]) > 0.0) {
911: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]));
912: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
913: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])));
914: } else {
915: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));
916: }
917: #else
918: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]));
919: #endif
921: /* U part */
922: for (j=a->diag[i+1]+1; j<a->diag[i]; j++) {
923: #if defined(PETSC_USE_COMPLEX)
924: if (PetscImaginaryPart(a->a[j]) > 0.0) {
925: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
926: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
927: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));
928: } else {
929: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
930: }
931: #else
932: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j]);
933: #endif
934: }
935: PetscViewerASCIIPrintf(viewer,"\n");
936: }
937: } else {
938: for (i=0; i<m; i++) {
939: PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i);
940: for (j=a->i[i]; j<a->i[i+1]; j++) {
941: #if defined(PETSC_USE_COMPLEX)
942: if (PetscImaginaryPart(a->a[j]) > 0.0) {
943: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
944: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
945: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));
946: } else {
947: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
948: }
949: #else
950: PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j]);
951: #endif
952: }
953: PetscViewerASCIIPrintf(viewer,"\n");
954: }
955: }
956: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
957: }
958: PetscViewerFlush(viewer);
959: return 0;
960: }
962: #include <petscdraw.h>
963: PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
964: {
965: Mat A = (Mat) Aa;
966: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
967: PetscInt i,j,m = A->rmap->n;
968: int color;
969: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
970: PetscViewer viewer;
971: PetscViewerFormat format;
972: const PetscScalar *aa;
974: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
975: PetscViewerGetFormat(viewer,&format);
976: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
978: /* loop over matrix elements drawing boxes */
979: MatSeqAIJGetArrayRead(A,&aa);
980: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
982: PetscDrawCollectiveBegin(draw);
983: /* Blue for negative, Cyan for zero and Red for positive */
984: color = PETSC_DRAW_BLUE;
985: for (i=0; i<m; i++) {
986: y_l = m - i - 1.0; y_r = y_l + 1.0;
987: for (j=a->i[i]; j<a->i[i+1]; j++) {
988: x_l = a->j[j]; x_r = x_l + 1.0;
989: if (PetscRealPart(aa[j]) >= 0.) continue;
990: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
991: }
992: }
993: color = PETSC_DRAW_CYAN;
994: for (i=0; i<m; i++) {
995: y_l = m - i - 1.0; y_r = y_l + 1.0;
996: for (j=a->i[i]; j<a->i[i+1]; j++) {
997: x_l = a->j[j]; x_r = x_l + 1.0;
998: if (aa[j] != 0.) continue;
999: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1000: }
1001: }
1002: color = PETSC_DRAW_RED;
1003: for (i=0; i<m; i++) {
1004: y_l = m - i - 1.0; y_r = y_l + 1.0;
1005: for (j=a->i[i]; j<a->i[i+1]; j++) {
1006: x_l = a->j[j]; x_r = x_l + 1.0;
1007: if (PetscRealPart(aa[j]) <= 0.) continue;
1008: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1009: }
1010: }
1011: PetscDrawCollectiveEnd(draw);
1012: } else {
1013: /* use contour shading to indicate magnitude of values */
1014: /* first determine max of all nonzero values */
1015: PetscReal minv = 0.0, maxv = 0.0;
1016: PetscInt nz = a->nz, count = 0;
1017: PetscDraw popup;
1020: for (i=0; i<nz; i++) {
1021: if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]);
1022: }
1023: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1024: PetscDrawGetPopup(draw,&popup);
1025: PetscDrawScalePopup(popup,minv,maxv);
1027: PetscDrawCollectiveBegin(draw);
1028: for (i=0; i<m; i++) {
1029: y_l = m - i - 1.0;
1030: y_r = y_l + 1.0;
1031: for (j=a->i[i]; j<a->i[i+1]; j++) {
1032: x_l = a->j[j];
1033: x_r = x_l + 1.0;
1034: color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv);
1035: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1036: count++;
1037: }
1038: }
1039: PetscDrawCollectiveEnd(draw);
1040: }
1041: MatSeqAIJRestoreArrayRead(A,&aa);
1042: return 0;
1043: }
1045: #include <petscdraw.h>
1046: PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer)
1047: {
1048: PetscDraw draw;
1049: PetscReal xr,yr,xl,yl,h,w;
1050: PetscBool isnull;
1052: PetscViewerDrawGetDraw(viewer,0,&draw);
1053: PetscDrawIsNull(draw,&isnull);
1054: if (isnull) return 0;
1056: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1057: xr += w; yr += h; xl = -w; yl = -h;
1058: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1059: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1060: PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);
1061: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1062: PetscDrawSave(draw);
1063: return 0;
1064: }
1066: PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer)
1067: {
1068: PetscBool iascii,isbinary,isdraw;
1070: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1071: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1072: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1073: if (iascii) {
1074: MatView_SeqAIJ_ASCII(A,viewer);
1075: } else if (isbinary) {
1076: MatView_SeqAIJ_Binary(A,viewer);
1077: } else if (isdraw) {
1078: MatView_SeqAIJ_Draw(A,viewer);
1079: }
1080: MatView_SeqAIJ_Inode(A,viewer);
1081: return 0;
1082: }
1084: PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode)
1085: {
1086: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1087: PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
1088: PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0;
1089: MatScalar *aa = a->a,*ap;
1090: PetscReal ratio = 0.6;
1092: if (mode == MAT_FLUSH_ASSEMBLY) return 0;
1093: MatSeqAIJInvalidateDiagonal(A);
1094: if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) {
1095: /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */
1096: MatAssemblyEnd_SeqAIJ_Inode(A,mode);
1097: return 0;
1098: }
1100: if (m) rmax = ailen[0]; /* determine row with most nonzeros */
1101: for (i=1; i<m; i++) {
1102: /* move each row back by the amount of empty slots (fshift) before it*/
1103: fshift += imax[i-1] - ailen[i-1];
1104: rmax = PetscMax(rmax,ailen[i]);
1105: if (fshift) {
1106: ip = aj + ai[i];
1107: ap = aa + ai[i];
1108: N = ailen[i];
1109: PetscArraymove(ip-fshift,ip,N);
1110: if (!A->structure_only) {
1111: PetscArraymove(ap-fshift,ap,N);
1112: }
1113: }
1114: ai[i] = ai[i-1] + ailen[i-1];
1115: }
1116: if (m) {
1117: fshift += imax[m-1] - ailen[m-1];
1118: ai[m] = ai[m-1] + ailen[m-1];
1119: }
1121: /* reset ilen and imax for each row */
1122: a->nonzerorowcnt = 0;
1123: if (A->structure_only) {
1124: PetscFree(a->imax);
1125: PetscFree(a->ilen);
1126: } else { /* !A->structure_only */
1127: for (i=0; i<m; i++) {
1128: ailen[i] = imax[i] = ai[i+1] - ai[i];
1129: a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0);
1130: }
1131: }
1132: a->nz = ai[m];
1135: MatMarkDiagonal_SeqAIJ(A);
1136: PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz);
1137: PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs);
1138: PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax);
1140: A->info.mallocs += a->reallocs;
1141: a->reallocs = 0;
1142: A->info.nz_unneeded = (PetscReal)fshift;
1143: a->rmax = rmax;
1145: if (!A->structure_only) {
1146: MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);
1147: }
1148: MatAssemblyEnd_SeqAIJ_Inode(A,mode);
1149: return 0;
1150: }
1152: PetscErrorCode MatRealPart_SeqAIJ(Mat A)
1153: {
1154: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1155: PetscInt i,nz = a->nz;
1156: MatScalar *aa;
1158: MatSeqAIJGetArray(A,&aa);
1159: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1160: MatSeqAIJRestoreArray(A,&aa);
1161: MatSeqAIJInvalidateDiagonal(A);
1162: return 0;
1163: }
1165: PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A)
1166: {
1167: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1168: PetscInt i,nz = a->nz;
1169: MatScalar *aa;
1171: MatSeqAIJGetArray(A,&aa);
1172: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1173: MatSeqAIJRestoreArray(A,&aa);
1174: MatSeqAIJInvalidateDiagonal(A);
1175: return 0;
1176: }
1178: PetscErrorCode MatZeroEntries_SeqAIJ(Mat A)
1179: {
1180: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1181: MatScalar *aa;
1183: MatSeqAIJGetArrayWrite(A,&aa);
1184: PetscArrayzero(aa,a->i[A->rmap->n]);
1185: MatSeqAIJRestoreArrayWrite(A,&aa);
1186: MatSeqAIJInvalidateDiagonal(A);
1187: return 0;
1188: }
1190: PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A)
1191: {
1192: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1194: PetscFree(a->perm);
1195: PetscFree(a->jmap);
1196: return 0;
1197: }
1199: PetscErrorCode MatDestroy_SeqAIJ(Mat A)
1200: {
1201: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1203: #if defined(PETSC_USE_LOG)
1204: PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz);
1205: #endif
1206: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
1207: MatResetPreallocationCOO_SeqAIJ(A);
1208: ISDestroy(&a->row);
1209: ISDestroy(&a->col);
1210: PetscFree(a->diag);
1211: PetscFree(a->ibdiag);
1212: PetscFree(a->imax);
1213: PetscFree(a->ilen);
1214: PetscFree(a->ipre);
1215: PetscFree3(a->idiag,a->mdiag,a->ssor_work);
1216: PetscFree(a->solve_work);
1217: ISDestroy(&a->icol);
1218: PetscFree(a->saved_values);
1219: PetscFree2(a->compressedrow.i,a->compressedrow.rindex);
1220: MatDestroy_SeqAIJ_Inode(A);
1221: PetscFree(A->data);
1223: /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this.
1224: That function is so heavily used (sometimes in an hidden way through multnumeric function pointers)
1225: that is hard to properly add this data to the MatProduct data. We free it here to avoid
1226: users reusing the matrix object with different data to incur in obscure segmentation faults
1227: due to different matrix sizes */
1228: PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL);
1230: PetscObjectChangeTypeName((PetscObject)A,NULL);
1231: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);
1232: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
1233: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
1234: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);
1235: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);
1236: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);
1237: #if defined(PETSC_HAVE_CUDA)
1238: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL);
1239: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL);
1240: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL);
1241: #endif
1242: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
1243: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL);
1244: #endif
1245: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL);
1246: #if defined(PETSC_HAVE_ELEMENTAL)
1247: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL);
1248: #endif
1249: #if defined(PETSC_HAVE_SCALAPACK)
1250: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL);
1251: #endif
1252: #if defined(PETSC_HAVE_HYPRE)
1253: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL);
1254: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL);
1255: #endif
1256: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL);
1257: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL);
1258: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL);
1259: PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);
1260: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);
1261: PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL);
1262: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);
1263: PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);
1264: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL);
1265: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL);
1266: PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL);
1267: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL);
1268: PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL);
1269: PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL);
1270: return 0;
1271: }
1273: PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg)
1274: {
1275: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1277: switch (op) {
1278: case MAT_ROW_ORIENTED:
1279: a->roworiented = flg;
1280: break;
1281: case MAT_KEEP_NONZERO_PATTERN:
1282: a->keepnonzeropattern = flg;
1283: break;
1284: case MAT_NEW_NONZERO_LOCATIONS:
1285: a->nonew = (flg ? 0 : 1);
1286: break;
1287: case MAT_NEW_NONZERO_LOCATION_ERR:
1288: a->nonew = (flg ? -1 : 0);
1289: break;
1290: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1291: a->nonew = (flg ? -2 : 0);
1292: break;
1293: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1294: a->nounused = (flg ? -1 : 0);
1295: break;
1296: case MAT_IGNORE_ZERO_ENTRIES:
1297: a->ignorezeroentries = flg;
1298: break;
1299: case MAT_SPD:
1300: case MAT_SYMMETRIC:
1301: case MAT_STRUCTURALLY_SYMMETRIC:
1302: case MAT_HERMITIAN:
1303: case MAT_SYMMETRY_ETERNAL:
1304: case MAT_STRUCTURE_ONLY:
1305: /* These options are handled directly by MatSetOption() */
1306: break;
1307: case MAT_FORCE_DIAGONAL_ENTRIES:
1308: case MAT_IGNORE_OFF_PROC_ENTRIES:
1309: case MAT_USE_HASH_TABLE:
1310: PetscInfo(A,"Option %s ignored\n",MatOptions[op]);
1311: break;
1312: case MAT_USE_INODES:
1313: MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg);
1314: break;
1315: case MAT_SUBMAT_SINGLEIS:
1316: A->submat_singleis = flg;
1317: break;
1318: case MAT_SORTED_FULL:
1319: if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull;
1320: else A->ops->setvalues = MatSetValues_SeqAIJ;
1321: break;
1322: case MAT_FORM_EXPLICIT_TRANSPOSE:
1323: A->form_explicit_transpose = flg;
1324: break;
1325: default:
1326: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
1327: }
1328: return 0;
1329: }
1331: PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v)
1332: {
1333: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1334: PetscInt i,j,n,*ai=a->i,*aj=a->j;
1335: PetscScalar *x;
1336: const PetscScalar *aa;
1338: VecGetLocalSize(v,&n);
1340: MatSeqAIJGetArrayRead(A,&aa);
1341: if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) {
1342: PetscInt *diag=a->diag;
1343: VecGetArrayWrite(v,&x);
1344: for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]];
1345: VecRestoreArrayWrite(v,&x);
1346: MatSeqAIJRestoreArrayRead(A,&aa);
1347: return 0;
1348: }
1350: VecGetArrayWrite(v,&x);
1351: for (i=0; i<n; i++) {
1352: x[i] = 0.0;
1353: for (j=ai[i]; j<ai[i+1]; j++) {
1354: if (aj[j] == i) {
1355: x[i] = aa[j];
1356: break;
1357: }
1358: }
1359: }
1360: VecRestoreArrayWrite(v,&x);
1361: MatSeqAIJRestoreArrayRead(A,&aa);
1362: return 0;
1363: }
1365: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1366: PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy)
1367: {
1368: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1369: const MatScalar *aa;
1370: PetscScalar *y;
1371: const PetscScalar *x;
1372: PetscInt m = A->rmap->n;
1373: #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1374: const MatScalar *v;
1375: PetscScalar alpha;
1376: PetscInt n,i,j;
1377: const PetscInt *idx,*ii,*ridx=NULL;
1378: Mat_CompressedRow cprow = a->compressedrow;
1379: PetscBool usecprow = cprow.use;
1380: #endif
1382: if (zz != yy) VecCopy(zz,yy);
1383: VecGetArrayRead(xx,&x);
1384: VecGetArray(yy,&y);
1385: MatSeqAIJGetArrayRead(A,&aa);
1387: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1388: fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y);
1389: #else
1390: if (usecprow) {
1391: m = cprow.nrows;
1392: ii = cprow.i;
1393: ridx = cprow.rindex;
1394: } else {
1395: ii = a->i;
1396: }
1397: for (i=0; i<m; i++) {
1398: idx = a->j + ii[i];
1399: v = aa + ii[i];
1400: n = ii[i+1] - ii[i];
1401: if (usecprow) {
1402: alpha = x[ridx[i]];
1403: } else {
1404: alpha = x[i];
1405: }
1406: for (j=0; j<n; j++) y[idx[j]] += alpha*v[j];
1407: }
1408: #endif
1409: PetscLogFlops(2.0*a->nz);
1410: VecRestoreArrayRead(xx,&x);
1411: VecRestoreArray(yy,&y);
1412: MatSeqAIJRestoreArrayRead(A,&aa);
1413: return 0;
1414: }
1416: PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy)
1417: {
1418: VecSet(yy,0.0);
1419: MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);
1420: return 0;
1421: }
1423: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1425: PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy)
1426: {
1427: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1428: PetscScalar *y;
1429: const PetscScalar *x;
1430: const MatScalar *aa,*a_a;
1431: PetscInt m=A->rmap->n;
1432: const PetscInt *aj,*ii,*ridx=NULL;
1433: PetscInt n,i;
1434: PetscScalar sum;
1435: PetscBool usecprow=a->compressedrow.use;
1437: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1438: #pragma disjoint(*x,*y,*aa)
1439: #endif
1441: if (a->inode.use && a->inode.checked) {
1442: MatMult_SeqAIJ_Inode(A,xx,yy);
1443: return 0;
1444: }
1445: MatSeqAIJGetArrayRead(A,&a_a);
1446: VecGetArrayRead(xx,&x);
1447: VecGetArray(yy,&y);
1448: ii = a->i;
1449: if (usecprow) { /* use compressed row format */
1450: PetscArrayzero(y,m);
1451: m = a->compressedrow.nrows;
1452: ii = a->compressedrow.i;
1453: ridx = a->compressedrow.rindex;
1454: for (i=0; i<m; i++) {
1455: n = ii[i+1] - ii[i];
1456: aj = a->j + ii[i];
1457: aa = a_a + ii[i];
1458: sum = 0.0;
1459: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1460: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1461: y[*ridx++] = sum;
1462: }
1463: } else { /* do not use compressed row format */
1464: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ)
1465: aj = a->j;
1466: aa = a_a;
1467: fortranmultaij_(&m,x,ii,aj,aa,y);
1468: #else
1469: for (i=0; i<m; i++) {
1470: n = ii[i+1] - ii[i];
1471: aj = a->j + ii[i];
1472: aa = a_a + ii[i];
1473: sum = 0.0;
1474: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1475: y[i] = sum;
1476: }
1477: #endif
1478: }
1479: PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);
1480: VecRestoreArrayRead(xx,&x);
1481: VecRestoreArray(yy,&y);
1482: MatSeqAIJRestoreArrayRead(A,&a_a);
1483: return 0;
1484: }
1486: PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy)
1487: {
1488: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1489: PetscScalar *y;
1490: const PetscScalar *x;
1491: const MatScalar *aa,*a_a;
1492: PetscInt m=A->rmap->n;
1493: const PetscInt *aj,*ii,*ridx=NULL;
1494: PetscInt n,i,nonzerorow=0;
1495: PetscScalar sum;
1496: PetscBool usecprow=a->compressedrow.use;
1498: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1499: #pragma disjoint(*x,*y,*aa)
1500: #endif
1502: MatSeqAIJGetArrayRead(A,&a_a);
1503: VecGetArrayRead(xx,&x);
1504: VecGetArray(yy,&y);
1505: if (usecprow) { /* use compressed row format */
1506: m = a->compressedrow.nrows;
1507: ii = a->compressedrow.i;
1508: ridx = a->compressedrow.rindex;
1509: for (i=0; i<m; i++) {
1510: n = ii[i+1] - ii[i];
1511: aj = a->j + ii[i];
1512: aa = a_a + ii[i];
1513: sum = 0.0;
1514: nonzerorow += (n>0);
1515: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1516: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1517: y[*ridx++] = sum;
1518: }
1519: } else { /* do not use compressed row format */
1520: ii = a->i;
1521: for (i=0; i<m; i++) {
1522: n = ii[i+1] - ii[i];
1523: aj = a->j + ii[i];
1524: aa = a_a + ii[i];
1525: sum = 0.0;
1526: nonzerorow += (n>0);
1527: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1528: y[i] = sum;
1529: }
1530: }
1531: PetscLogFlops(2.0*a->nz - nonzerorow);
1532: VecRestoreArrayRead(xx,&x);
1533: VecRestoreArray(yy,&y);
1534: MatSeqAIJRestoreArrayRead(A,&a_a);
1535: return 0;
1536: }
1538: PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz)
1539: {
1540: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1541: PetscScalar *y,*z;
1542: const PetscScalar *x;
1543: const MatScalar *aa,*a_a;
1544: PetscInt m = A->rmap->n,*aj,*ii;
1545: PetscInt n,i,*ridx=NULL;
1546: PetscScalar sum;
1547: PetscBool usecprow=a->compressedrow.use;
1549: MatSeqAIJGetArrayRead(A,&a_a);
1550: VecGetArrayRead(xx,&x);
1551: VecGetArrayPair(yy,zz,&y,&z);
1552: if (usecprow) { /* use compressed row format */
1553: if (zz != yy) {
1554: PetscArraycpy(z,y,m);
1555: }
1556: m = a->compressedrow.nrows;
1557: ii = a->compressedrow.i;
1558: ridx = a->compressedrow.rindex;
1559: for (i=0; i<m; i++) {
1560: n = ii[i+1] - ii[i];
1561: aj = a->j + ii[i];
1562: aa = a_a + ii[i];
1563: sum = y[*ridx];
1564: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1565: z[*ridx++] = sum;
1566: }
1567: } else { /* do not use compressed row format */
1568: ii = a->i;
1569: for (i=0; i<m; i++) {
1570: n = ii[i+1] - ii[i];
1571: aj = a->j + ii[i];
1572: aa = a_a + ii[i];
1573: sum = y[i];
1574: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1575: z[i] = sum;
1576: }
1577: }
1578: PetscLogFlops(2.0*a->nz);
1579: VecRestoreArrayRead(xx,&x);
1580: VecRestoreArrayPair(yy,zz,&y,&z);
1581: MatSeqAIJRestoreArrayRead(A,&a_a);
1582: return 0;
1583: }
1585: #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h>
1586: PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz)
1587: {
1588: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1589: PetscScalar *y,*z;
1590: const PetscScalar *x;
1591: const MatScalar *aa,*a_a;
1592: const PetscInt *aj,*ii,*ridx=NULL;
1593: PetscInt m = A->rmap->n,n,i;
1594: PetscScalar sum;
1595: PetscBool usecprow=a->compressedrow.use;
1597: if (a->inode.use && a->inode.checked) {
1598: MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz);
1599: return 0;
1600: }
1601: MatSeqAIJGetArrayRead(A,&a_a);
1602: VecGetArrayRead(xx,&x);
1603: VecGetArrayPair(yy,zz,&y,&z);
1604: if (usecprow) { /* use compressed row format */
1605: if (zz != yy) {
1606: PetscArraycpy(z,y,m);
1607: }
1608: m = a->compressedrow.nrows;
1609: ii = a->compressedrow.i;
1610: ridx = a->compressedrow.rindex;
1611: for (i=0; i<m; i++) {
1612: n = ii[i+1] - ii[i];
1613: aj = a->j + ii[i];
1614: aa = a_a + ii[i];
1615: sum = y[*ridx];
1616: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1617: z[*ridx++] = sum;
1618: }
1619: } else { /* do not use compressed row format */
1620: ii = a->i;
1621: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ)
1622: aj = a->j;
1623: aa = a_a;
1624: fortranmultaddaij_(&m,x,ii,aj,aa,y,z);
1625: #else
1626: for (i=0; i<m; i++) {
1627: n = ii[i+1] - ii[i];
1628: aj = a->j + ii[i];
1629: aa = a_a + ii[i];
1630: sum = y[i];
1631: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1632: z[i] = sum;
1633: }
1634: #endif
1635: }
1636: PetscLogFlops(2.0*a->nz);
1637: VecRestoreArrayRead(xx,&x);
1638: VecRestoreArrayPair(yy,zz,&y,&z);
1639: MatSeqAIJRestoreArrayRead(A,&a_a);
1640: return 0;
1641: }
1643: /*
1644: Adds diagonal pointers to sparse matrix structure.
1645: */
1646: PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A)
1647: {
1648: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1649: PetscInt i,j,m = A->rmap->n;
1650: PetscBool alreadySet = PETSC_TRUE;
1652: if (!a->diag) {
1653: PetscMalloc1(m,&a->diag);
1654: PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));
1655: alreadySet = PETSC_FALSE;
1656: }
1657: for (i=0; i<A->rmap->n; i++) {
1658: /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */
1659: if (alreadySet) {
1660: PetscInt pos = a->diag[i];
1661: if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue;
1662: }
1664: a->diag[i] = a->i[i+1];
1665: for (j=a->i[i]; j<a->i[i+1]; j++) {
1666: if (a->j[j] == i) {
1667: a->diag[i] = j;
1668: break;
1669: }
1670: }
1671: }
1672: return 0;
1673: }
1675: PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v)
1676: {
1677: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1678: const PetscInt *diag = (const PetscInt*)a->diag;
1679: const PetscInt *ii = (const PetscInt*) a->i;
1680: PetscInt i,*mdiag = NULL;
1681: PetscInt cnt = 0; /* how many diagonals are missing */
1683: if (!A->preallocated || !a->nz) {
1684: MatSeqAIJSetPreallocation(A,1,NULL);
1685: MatShift_Basic(A,v);
1686: return 0;
1687: }
1689: if (a->diagonaldense) {
1690: cnt = 0;
1691: } else {
1692: PetscCalloc1(A->rmap->n,&mdiag);
1693: for (i=0; i<A->rmap->n; i++) {
1694: if (diag[i] >= ii[i+1]) {
1695: cnt++;
1696: mdiag[i] = 1;
1697: }
1698: }
1699: }
1700: if (!cnt) {
1701: MatShift_Basic(A,v);
1702: } else {
1703: PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */
1704: PetscInt *oldj = a->j, *oldi = a->i;
1705: PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij;
1707: a->a = NULL;
1708: a->j = NULL;
1709: a->i = NULL;
1710: /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */
1711: for (i=0; i<A->rmap->n; i++) {
1712: a->imax[i] += mdiag[i];
1713: a->imax[i] = PetscMin(a->imax[i],A->cmap->n);
1714: }
1715: MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax);
1717: /* copy old values into new matrix data structure */
1718: for (i=0; i<A->rmap->n; i++) {
1719: MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES);
1720: if (i < A->cmap->n) {
1721: MatSetValue(A,i,i,v,ADD_VALUES);
1722: }
1723: }
1724: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1725: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1726: if (singlemalloc) {
1727: PetscFree3(olda,oldj,oldi);
1728: } else {
1729: if (free_a) PetscFree(olda);
1730: if (free_ij) PetscFree(oldj);
1731: if (free_ij) PetscFree(oldi);
1732: }
1733: }
1734: PetscFree(mdiag);
1735: a->diagonaldense = PETSC_TRUE;
1736: return 0;
1737: }
1739: /*
1740: Checks for missing diagonals
1741: */
1742: PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d)
1743: {
1744: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1745: PetscInt *diag,*ii = a->i,i;
1747: *missing = PETSC_FALSE;
1748: if (A->rmap->n > 0 && !ii) {
1749: *missing = PETSC_TRUE;
1750: if (d) *d = 0;
1751: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
1752: } else {
1753: PetscInt n;
1754: n = PetscMin(A->rmap->n, A->cmap->n);
1755: diag = a->diag;
1756: for (i=0; i<n; i++) {
1757: if (diag[i] >= ii[i+1]) {
1758: *missing = PETSC_TRUE;
1759: if (d) *d = i;
1760: PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i);
1761: break;
1762: }
1763: }
1764: }
1765: return 0;
1766: }
1768: #include <petscblaslapack.h>
1769: #include <petsc/private/kernels/blockinvert.h>
1771: /*
1772: Note that values is allocated externally by the PC and then passed into this routine
1773: */
1774: PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag)
1775: {
1776: PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots;
1777: PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
1778: const PetscReal shift = 0.0;
1779: PetscInt ipvt[5];
1780: PetscScalar work[25],*v_work;
1782: allowzeropivot = PetscNot(A->erroriffailure);
1783: for (i=0; i<nblocks; i++) ncnt += bsizes[i];
1785: for (i=0; i<nblocks; i++) {
1786: bsizemax = PetscMax(bsizemax,bsizes[i]);
1787: }
1788: PetscMalloc1(bsizemax,&indx);
1789: if (bsizemax > 7) {
1790: PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots);
1791: }
1792: ncnt = 0;
1793: for (i=0; i<nblocks; i++) {
1794: for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j;
1795: MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag);
1796: switch (bsizes[i]) {
1797: case 1:
1798: *diag = 1.0/(*diag);
1799: break;
1800: case 2:
1801: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
1802: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1803: PetscKernel_A_gets_transpose_A_2(diag);
1804: break;
1805: case 3:
1806: PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
1807: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1808: PetscKernel_A_gets_transpose_A_3(diag);
1809: break;
1810: case 4:
1811: PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
1812: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1813: PetscKernel_A_gets_transpose_A_4(diag);
1814: break;
1815: case 5:
1816: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
1817: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1818: PetscKernel_A_gets_transpose_A_5(diag);
1819: break;
1820: case 6:
1821: PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
1822: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1823: PetscKernel_A_gets_transpose_A_6(diag);
1824: break;
1825: case 7:
1826: PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
1827: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1828: PetscKernel_A_gets_transpose_A_7(diag);
1829: break;
1830: default:
1831: PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
1832: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1833: PetscKernel_A_gets_transpose_A_N(diag,bsizes[i]);
1834: }
1835: ncnt += bsizes[i];
1836: diag += bsizes[i]*bsizes[i];
1837: }
1838: if (bsizemax > 7) {
1839: PetscFree2(v_work,v_pivots);
1840: }
1841: PetscFree(indx);
1842: return 0;
1843: }
1845: /*
1846: Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways
1847: */
1848: PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift)
1849: {
1850: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
1851: PetscInt i,*diag,m = A->rmap->n;
1852: const MatScalar *v;
1853: PetscScalar *idiag,*mdiag;
1855: if (a->idiagvalid) return 0;
1856: MatMarkDiagonal_SeqAIJ(A);
1857: diag = a->diag;
1858: if (!a->idiag) {
1859: PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);
1860: PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar));
1861: }
1863: mdiag = a->mdiag;
1864: idiag = a->idiag;
1865: MatSeqAIJGetArrayRead(A,&v);
1866: if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) {
1867: for (i=0; i<m; i++) {
1868: mdiag[i] = v[diag[i]];
1869: if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */
1870: if (PetscRealPart(fshift)) {
1871: PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i);
1872: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1873: A->factorerror_zeropivot_value = 0.0;
1874: A->factorerror_zeropivot_row = i;
1875: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i);
1876: }
1877: idiag[i] = 1.0/v[diag[i]];
1878: }
1879: PetscLogFlops(m);
1880: } else {
1881: for (i=0; i<m; i++) {
1882: mdiag[i] = v[diag[i]];
1883: idiag[i] = omega/(fshift + v[diag[i]]);
1884: }
1885: PetscLogFlops(2.0*m);
1886: }
1887: a->idiagvalid = PETSC_TRUE;
1888: MatSeqAIJRestoreArrayRead(A,&v);
1889: return 0;
1890: }
1892: #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h>
1893: PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1894: {
1895: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1896: PetscScalar *x,d,sum,*t,scale;
1897: const MatScalar *v,*idiag=NULL,*mdiag,*aa;
1898: const PetscScalar *b, *bs,*xb, *ts;
1899: PetscInt n,m = A->rmap->n,i;
1900: const PetscInt *idx,*diag;
1902: if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) {
1903: MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx);
1904: return 0;
1905: }
1906: its = its*lits;
1908: if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */
1909: if (!a->idiagvalid) MatInvertDiagonal_SeqAIJ(A,omega,fshift);
1910: a->fshift = fshift;
1911: a->omega = omega;
1913: diag = a->diag;
1914: t = a->ssor_work;
1915: idiag = a->idiag;
1916: mdiag = a->mdiag;
1918: MatSeqAIJGetArrayRead(A,&aa);
1919: VecGetArray(xx,&x);
1920: VecGetArrayRead(bb,&b);
1921: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
1922: if (flag == SOR_APPLY_UPPER) {
1923: /* apply (U + D/omega) to the vector */
1924: bs = b;
1925: for (i=0; i<m; i++) {
1926: d = fshift + mdiag[i];
1927: n = a->i[i+1] - diag[i] - 1;
1928: idx = a->j + diag[i] + 1;
1929: v = aa + diag[i] + 1;
1930: sum = b[i]*d/omega;
1931: PetscSparseDensePlusDot(sum,bs,v,idx,n);
1932: x[i] = sum;
1933: }
1934: VecRestoreArray(xx,&x);
1935: VecRestoreArrayRead(bb,&b);
1936: MatSeqAIJRestoreArrayRead(A,&aa);
1937: PetscLogFlops(a->nz);
1938: return 0;
1939: }
1942: else if (flag & SOR_EISENSTAT) {
1943: /* Let A = L + U + D; where L is lower triangular,
1944: U is upper triangular, E = D/omega; This routine applies
1946: (L + E)^{-1} A (U + E)^{-1}
1948: to a vector efficiently using Eisenstat's trick.
1949: */
1950: scale = (2.0/omega) - 1.0;
1952: /* x = (E + U)^{-1} b */
1953: for (i=m-1; i>=0; i--) {
1954: n = a->i[i+1] - diag[i] - 1;
1955: idx = a->j + diag[i] + 1;
1956: v = aa + diag[i] + 1;
1957: sum = b[i];
1958: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1959: x[i] = sum*idiag[i];
1960: }
1962: /* t = b - (2*E - D)x */
1963: v = aa;
1964: for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i];
1966: /* t = (E + L)^{-1}t */
1967: ts = t;
1968: diag = a->diag;
1969: for (i=0; i<m; i++) {
1970: n = diag[i] - a->i[i];
1971: idx = a->j + a->i[i];
1972: v = aa + a->i[i];
1973: sum = t[i];
1974: PetscSparseDenseMinusDot(sum,ts,v,idx,n);
1975: t[i] = sum*idiag[i];
1976: /* x = x + t */
1977: x[i] += t[i];
1978: }
1980: PetscLogFlops(6.0*m-1 + 2.0*a->nz);
1981: VecRestoreArray(xx,&x);
1982: VecRestoreArrayRead(bb,&b);
1983: return 0;
1984: }
1985: if (flag & SOR_ZERO_INITIAL_GUESS) {
1986: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1987: for (i=0; i<m; i++) {
1988: n = diag[i] - a->i[i];
1989: idx = a->j + a->i[i];
1990: v = aa + a->i[i];
1991: sum = b[i];
1992: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1993: t[i] = sum;
1994: x[i] = sum*idiag[i];
1995: }
1996: xb = t;
1997: PetscLogFlops(a->nz);
1998: } else xb = b;
1999: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
2000: for (i=m-1; i>=0; i--) {
2001: n = a->i[i+1] - diag[i] - 1;
2002: idx = a->j + diag[i] + 1;
2003: v = aa + diag[i] + 1;
2004: sum = xb[i];
2005: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2006: if (xb == b) {
2007: x[i] = sum*idiag[i];
2008: } else {
2009: x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2010: }
2011: }
2012: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
2013: }
2014: its--;
2015: }
2016: while (its--) {
2017: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
2018: for (i=0; i<m; i++) {
2019: /* lower */
2020: n = diag[i] - a->i[i];
2021: idx = a->j + a->i[i];
2022: v = aa + a->i[i];
2023: sum = b[i];
2024: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2025: t[i] = sum; /* save application of the lower-triangular part */
2026: /* upper */
2027: n = a->i[i+1] - diag[i] - 1;
2028: idx = a->j + diag[i] + 1;
2029: v = aa + diag[i] + 1;
2030: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2031: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2032: }
2033: xb = t;
2034: PetscLogFlops(2.0*a->nz);
2035: } else xb = b;
2036: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
2037: for (i=m-1; i>=0; i--) {
2038: sum = xb[i];
2039: if (xb == b) {
2040: /* whole matrix (no checkpointing available) */
2041: n = a->i[i+1] - a->i[i];
2042: idx = a->j + a->i[i];
2043: v = aa + a->i[i];
2044: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2045: x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i];
2046: } else { /* lower-triangular part has been saved, so only apply upper-triangular */
2047: n = a->i[i+1] - diag[i] - 1;
2048: idx = a->j + diag[i] + 1;
2049: v = aa + diag[i] + 1;
2050: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2051: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2052: }
2053: }
2054: if (xb == b) {
2055: PetscLogFlops(2.0*a->nz);
2056: } else {
2057: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
2058: }
2059: }
2060: }
2061: MatSeqAIJRestoreArrayRead(A,&aa);
2062: VecRestoreArray(xx,&x);
2063: VecRestoreArrayRead(bb,&b);
2064: return 0;
2065: }
2067: PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info)
2068: {
2069: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2071: info->block_size = 1.0;
2072: info->nz_allocated = a->maxnz;
2073: info->nz_used = a->nz;
2074: info->nz_unneeded = (a->maxnz - a->nz);
2075: info->assemblies = A->num_ass;
2076: info->mallocs = A->info.mallocs;
2077: info->memory = ((PetscObject)A)->mem;
2078: if (A->factortype) {
2079: info->fill_ratio_given = A->info.fill_ratio_given;
2080: info->fill_ratio_needed = A->info.fill_ratio_needed;
2081: info->factor_mallocs = A->info.factor_mallocs;
2082: } else {
2083: info->fill_ratio_given = 0;
2084: info->fill_ratio_needed = 0;
2085: info->factor_mallocs = 0;
2086: }
2087: return 0;
2088: }
2090: PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2091: {
2092: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2093: PetscInt i,m = A->rmap->n - 1;
2094: const PetscScalar *xx;
2095: PetscScalar *bb,*aa;
2096: PetscInt d = 0;
2098: if (x && b) {
2099: VecGetArrayRead(x,&xx);
2100: VecGetArray(b,&bb);
2101: for (i=0; i<N; i++) {
2103: if (rows[i] >= A->cmap->n) continue;
2104: bb[rows[i]] = diag*xx[rows[i]];
2105: }
2106: VecRestoreArrayRead(x,&xx);
2107: VecRestoreArray(b,&bb);
2108: }
2110: MatSeqAIJGetArray(A,&aa);
2111: if (a->keepnonzeropattern) {
2112: for (i=0; i<N; i++) {
2114: PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]]);
2115: }
2116: if (diag != 0.0) {
2117: for (i=0; i<N; i++) {
2118: d = rows[i];
2119: if (rows[i] >= A->cmap->n) continue;
2121: }
2122: for (i=0; i<N; i++) {
2123: if (rows[i] >= A->cmap->n) continue;
2124: aa[a->diag[rows[i]]] = diag;
2125: }
2126: }
2127: } else {
2128: if (diag != 0.0) {
2129: for (i=0; i<N; i++) {
2131: if (a->ilen[rows[i]] > 0) {
2132: if (rows[i] >= A->cmap->n) {
2133: a->ilen[rows[i]] = 0;
2134: } else {
2135: a->ilen[rows[i]] = 1;
2136: aa[a->i[rows[i]]] = diag;
2137: a->j[a->i[rows[i]]] = rows[i];
2138: }
2139: } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */
2140: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
2141: }
2142: }
2143: } else {
2144: for (i=0; i<N; i++) {
2146: a->ilen[rows[i]] = 0;
2147: }
2148: }
2149: A->nonzerostate++;
2150: }
2151: MatSeqAIJRestoreArray(A,&aa);
2152: (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);
2153: return 0;
2154: }
2156: PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2157: {
2158: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2159: PetscInt i,j,m = A->rmap->n - 1,d = 0;
2160: PetscBool missing,*zeroed,vecs = PETSC_FALSE;
2161: const PetscScalar *xx;
2162: PetscScalar *bb,*aa;
2164: if (!N) return 0;
2165: MatSeqAIJGetArray(A,&aa);
2166: if (x && b) {
2167: VecGetArrayRead(x,&xx);
2168: VecGetArray(b,&bb);
2169: vecs = PETSC_TRUE;
2170: }
2171: PetscCalloc1(A->rmap->n,&zeroed);
2172: for (i=0; i<N; i++) {
2174: PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]]);
2176: zeroed[rows[i]] = PETSC_TRUE;
2177: }
2178: for (i=0; i<A->rmap->n; i++) {
2179: if (!zeroed[i]) {
2180: for (j=a->i[i]; j<a->i[i+1]; j++) {
2181: if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) {
2182: if (vecs) bb[i] -= aa[j]*xx[a->j[j]];
2183: aa[j] = 0.0;
2184: }
2185: }
2186: } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i];
2187: }
2188: if (x && b) {
2189: VecRestoreArrayRead(x,&xx);
2190: VecRestoreArray(b,&bb);
2191: }
2192: PetscFree(zeroed);
2193: if (diag != 0.0) {
2194: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
2195: if (missing) {
2196: for (i=0; i<N; i++) {
2197: if (rows[i] >= A->cmap->N) continue;
2199: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
2200: }
2201: } else {
2202: for (i=0; i<N; i++) {
2203: aa[a->diag[rows[i]]] = diag;
2204: }
2205: }
2206: }
2207: MatSeqAIJRestoreArray(A,&aa);
2208: (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);
2209: return 0;
2210: }
2212: PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
2213: {
2214: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2215: const PetscScalar *aa;
2216: PetscInt *itmp;
2218: MatSeqAIJGetArrayRead(A,&aa);
2219: *nz = a->i[row+1] - a->i[row];
2220: if (v) *v = (PetscScalar*)(aa + a->i[row]);
2221: if (idx) {
2222: itmp = a->j + a->i[row];
2223: if (*nz) *idx = itmp;
2224: else *idx = NULL;
2225: }
2226: MatSeqAIJRestoreArrayRead(A,&aa);
2227: return 0;
2228: }
2230: PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
2231: {
2232: if (nz) *nz = 0;
2233: if (idx) *idx = NULL;
2234: if (v) *v = NULL;
2235: return 0;
2236: }
2238: PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm)
2239: {
2240: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2241: const MatScalar *v;
2242: PetscReal sum = 0.0;
2243: PetscInt i,j;
2245: MatSeqAIJGetArrayRead(A,&v);
2246: if (type == NORM_FROBENIUS) {
2247: #if defined(PETSC_USE_REAL___FP16)
2248: PetscBLASInt one = 1,nz = a->nz;
2249: PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one));
2250: #else
2251: for (i=0; i<a->nz; i++) {
2252: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
2253: }
2254: *nrm = PetscSqrtReal(sum);
2255: #endif
2256: PetscLogFlops(2.0*a->nz);
2257: } else if (type == NORM_1) {
2258: PetscReal *tmp;
2259: PetscInt *jj = a->j;
2260: PetscCalloc1(A->cmap->n+1,&tmp);
2261: *nrm = 0.0;
2262: for (j=0; j<a->nz; j++) {
2263: tmp[*jj++] += PetscAbsScalar(*v); v++;
2264: }
2265: for (j=0; j<A->cmap->n; j++) {
2266: if (tmp[j] > *nrm) *nrm = tmp[j];
2267: }
2268: PetscFree(tmp);
2269: PetscLogFlops(PetscMax(a->nz-1,0));
2270: } else if (type == NORM_INFINITY) {
2271: *nrm = 0.0;
2272: for (j=0; j<A->rmap->n; j++) {
2273: const PetscScalar *v2 = v + a->i[j];
2274: sum = 0.0;
2275: for (i=0; i<a->i[j+1]-a->i[j]; i++) {
2276: sum += PetscAbsScalar(*v2); v2++;
2277: }
2278: if (sum > *nrm) *nrm = sum;
2279: }
2280: PetscLogFlops(PetscMax(a->nz-1,0));
2281: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm");
2282: MatSeqAIJRestoreArrayRead(A,&v);
2283: return 0;
2284: }
2286: /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */
2287: PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B)
2288: {
2289: PetscInt i,j,anzj;
2290: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b;
2291: PetscInt an=A->cmap->N,am=A->rmap->N;
2292: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
2294: /* Allocate space for symbolic transpose info and work array */
2295: PetscCalloc1(an+1,&ati);
2296: PetscMalloc1(ai[am],&atj);
2297: PetscMalloc1(an,&atfill);
2299: /* Walk through aj and count ## of non-zeros in each row of A^T. */
2300: /* Note: offset by 1 for fast conversion into csr format. */
2301: for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1;
2302: /* Form ati for csr format of A^T. */
2303: for (i=0;i<an;i++) ati[i+1] += ati[i];
2305: /* Copy ati into atfill so we have locations of the next free space in atj */
2306: PetscArraycpy(atfill,ati,an);
2308: /* Walk through A row-wise and mark nonzero entries of A^T. */
2309: for (i=0;i<am;i++) {
2310: anzj = ai[i+1] - ai[i];
2311: for (j=0;j<anzj;j++) {
2312: atj[atfill[*aj]] = i;
2313: atfill[*aj++] += 1;
2314: }
2315: }
2317: /* Clean up temporary space and complete requests. */
2318: PetscFree(atfill);
2319: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);
2320: MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
2321: MatSetType(*B,((PetscObject)A)->type_name);
2323: b = (Mat_SeqAIJ*)((*B)->data);
2324: b->free_a = PETSC_FALSE;
2325: b->free_ij = PETSC_TRUE;
2326: b->nonew = 0;
2327: return 0;
2328: }
2330: PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2331: {
2332: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data;
2333: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2334: const MatScalar *va,*vb;
2335: PetscInt ma,na,mb,nb, i;
2337: MatGetSize(A,&ma,&na);
2338: MatGetSize(B,&mb,&nb);
2339: if (ma!=nb || na!=mb) {
2340: *f = PETSC_FALSE;
2341: return 0;
2342: }
2343: MatSeqAIJGetArrayRead(A,&va);
2344: MatSeqAIJGetArrayRead(B,&vb);
2345: aii = aij->i; bii = bij->i;
2346: adx = aij->j; bdx = bij->j;
2347: PetscMalloc1(ma,&aptr);
2348: PetscMalloc1(mb,&bptr);
2349: for (i=0; i<ma; i++) aptr[i] = aii[i];
2350: for (i=0; i<mb; i++) bptr[i] = bii[i];
2352: *f = PETSC_TRUE;
2353: for (i=0; i<ma; i++) {
2354: while (aptr[i]<aii[i+1]) {
2355: PetscInt idc,idr;
2356: PetscScalar vc,vr;
2357: /* column/row index/value */
2358: idc = adx[aptr[i]];
2359: idr = bdx[bptr[idc]];
2360: vc = va[aptr[i]];
2361: vr = vb[bptr[idc]];
2362: if (i!=idr || PetscAbsScalar(vc-vr) > tol) {
2363: *f = PETSC_FALSE;
2364: goto done;
2365: } else {
2366: aptr[i]++;
2367: if (B || i!=idc) bptr[idc]++;
2368: }
2369: }
2370: }
2371: done:
2372: PetscFree(aptr);
2373: PetscFree(bptr);
2374: MatSeqAIJRestoreArrayRead(A,&va);
2375: MatSeqAIJRestoreArrayRead(B,&vb);
2376: return 0;
2377: }
2379: PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2380: {
2381: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data;
2382: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2383: MatScalar *va,*vb;
2384: PetscInt ma,na,mb,nb, i;
2386: MatGetSize(A,&ma,&na);
2387: MatGetSize(B,&mb,&nb);
2388: if (ma!=nb || na!=mb) {
2389: *f = PETSC_FALSE;
2390: return 0;
2391: }
2392: aii = aij->i; bii = bij->i;
2393: adx = aij->j; bdx = bij->j;
2394: va = aij->a; vb = bij->a;
2395: PetscMalloc1(ma,&aptr);
2396: PetscMalloc1(mb,&bptr);
2397: for (i=0; i<ma; i++) aptr[i] = aii[i];
2398: for (i=0; i<mb; i++) bptr[i] = bii[i];
2400: *f = PETSC_TRUE;
2401: for (i=0; i<ma; i++) {
2402: while (aptr[i]<aii[i+1]) {
2403: PetscInt idc,idr;
2404: PetscScalar vc,vr;
2405: /* column/row index/value */
2406: idc = adx[aptr[i]];
2407: idr = bdx[bptr[idc]];
2408: vc = va[aptr[i]];
2409: vr = vb[bptr[idc]];
2410: if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) {
2411: *f = PETSC_FALSE;
2412: goto done;
2413: } else {
2414: aptr[i]++;
2415: if (B || i!=idc) bptr[idc]++;
2416: }
2417: }
2418: }
2419: done:
2420: PetscFree(aptr);
2421: PetscFree(bptr);
2422: return 0;
2423: }
2425: PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2426: {
2427: MatIsTranspose_SeqAIJ(A,A,tol,f);
2428: return 0;
2429: }
2431: PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2432: {
2433: MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);
2434: return 0;
2435: }
2437: PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr)
2438: {
2439: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2440: const PetscScalar *l,*r;
2441: PetscScalar x;
2442: MatScalar *v;
2443: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz;
2444: const PetscInt *jj;
2446: if (ll) {
2447: /* The local size is used so that VecMPI can be passed to this routine
2448: by MatDiagonalScale_MPIAIJ */
2449: VecGetLocalSize(ll,&m);
2451: VecGetArrayRead(ll,&l);
2452: MatSeqAIJGetArray(A,&v);
2453: for (i=0; i<m; i++) {
2454: x = l[i];
2455: M = a->i[i+1] - a->i[i];
2456: for (j=0; j<M; j++) (*v++) *= x;
2457: }
2458: VecRestoreArrayRead(ll,&l);
2459: PetscLogFlops(nz);
2460: MatSeqAIJRestoreArray(A,&v);
2461: }
2462: if (rr) {
2463: VecGetLocalSize(rr,&n);
2465: VecGetArrayRead(rr,&r);
2466: MatSeqAIJGetArray(A,&v);
2467: jj = a->j;
2468: for (i=0; i<nz; i++) (*v++) *= r[*jj++];
2469: MatSeqAIJRestoreArray(A,&v);
2470: VecRestoreArrayRead(rr,&r);
2471: PetscLogFlops(nz);
2472: }
2473: MatSeqAIJInvalidateDiagonal(A);
2474: return 0;
2475: }
2477: PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B)
2478: {
2479: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c;
2480: PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens;
2481: PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi;
2482: const PetscInt *irow,*icol;
2483: const PetscScalar *aa;
2484: PetscInt nrows,ncols;
2485: PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen;
2486: MatScalar *a_new,*mat_a;
2487: Mat C;
2488: PetscBool stride;
2490: ISGetIndices(isrow,&irow);
2491: ISGetLocalSize(isrow,&nrows);
2492: ISGetLocalSize(iscol,&ncols);
2494: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
2495: if (stride) {
2496: ISStrideGetInfo(iscol,&first,&step);
2497: } else {
2498: first = 0;
2499: step = 0;
2500: }
2501: if (stride && step == 1) {
2502: /* special case of contiguous rows */
2503: PetscMalloc2(nrows,&lens,nrows,&starts);
2504: /* loop over new rows determining lens and starting points */
2505: for (i=0; i<nrows; i++) {
2506: kstart = ai[irow[i]];
2507: kend = kstart + ailen[irow[i]];
2508: starts[i] = kstart;
2509: for (k=kstart; k<kend; k++) {
2510: if (aj[k] >= first) {
2511: starts[i] = k;
2512: break;
2513: }
2514: }
2515: sum = 0;
2516: while (k < kend) {
2517: if (aj[k++] >= first+ncols) break;
2518: sum++;
2519: }
2520: lens[i] = sum;
2521: }
2522: /* create submatrix */
2523: if (scall == MAT_REUSE_MATRIX) {
2524: PetscInt n_cols,n_rows;
2525: MatGetSize(*B,&n_rows,&n_cols);
2527: MatZeroEntries(*B);
2528: C = *B;
2529: } else {
2530: PetscInt rbs,cbs;
2531: MatCreate(PetscObjectComm((PetscObject)A),&C);
2532: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2533: ISGetBlockSize(isrow,&rbs);
2534: ISGetBlockSize(iscol,&cbs);
2535: MatSetBlockSizes(C,rbs,cbs);
2536: MatSetType(C,((PetscObject)A)->type_name);
2537: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2538: }
2539: c = (Mat_SeqAIJ*)C->data;
2541: /* loop over rows inserting into submatrix */
2542: a_new = c->a;
2543: j_new = c->j;
2544: i_new = c->i;
2545: MatSeqAIJGetArrayRead(A,&aa);
2546: for (i=0; i<nrows; i++) {
2547: ii = starts[i];
2548: lensi = lens[i];
2549: for (k=0; k<lensi; k++) {
2550: *j_new++ = aj[ii+k] - first;
2551: }
2552: PetscArraycpy(a_new,aa + starts[i],lensi);
2553: a_new += lensi;
2554: i_new[i+1] = i_new[i] + lensi;
2555: c->ilen[i] = lensi;
2556: }
2557: MatSeqAIJRestoreArrayRead(A,&aa);
2558: PetscFree2(lens,starts);
2559: } else {
2560: ISGetIndices(iscol,&icol);
2561: PetscCalloc1(oldcols,&smap);
2562: PetscMalloc1(1+nrows,&lens);
2563: for (i=0; i<ncols; i++) {
2565: smap[icol[i]] = i+1;
2566: }
2568: /* determine lens of each row */
2569: for (i=0; i<nrows; i++) {
2570: kstart = ai[irow[i]];
2571: kend = kstart + a->ilen[irow[i]];
2572: lens[i] = 0;
2573: for (k=kstart; k<kend; k++) {
2574: if (smap[aj[k]]) {
2575: lens[i]++;
2576: }
2577: }
2578: }
2579: /* Create and fill new matrix */
2580: if (scall == MAT_REUSE_MATRIX) {
2581: PetscBool equal;
2583: c = (Mat_SeqAIJ*)((*B)->data);
2585: PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal);
2587: PetscArrayzero(c->ilen,(*B)->rmap->n);
2588: C = *B;
2589: } else {
2590: PetscInt rbs,cbs;
2591: MatCreate(PetscObjectComm((PetscObject)A),&C);
2592: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2593: ISGetBlockSize(isrow,&rbs);
2594: ISGetBlockSize(iscol,&cbs);
2595: MatSetBlockSizes(C,rbs,cbs);
2596: MatSetType(C,((PetscObject)A)->type_name);
2597: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2598: }
2599: MatSeqAIJGetArrayRead(A,&aa);
2600: c = (Mat_SeqAIJ*)(C->data);
2601: for (i=0; i<nrows; i++) {
2602: row = irow[i];
2603: kstart = ai[row];
2604: kend = kstart + a->ilen[row];
2605: mat_i = c->i[i];
2606: mat_j = c->j + mat_i;
2607: mat_a = c->a + mat_i;
2608: mat_ilen = c->ilen + i;
2609: for (k=kstart; k<kend; k++) {
2610: if ((tcol=smap[a->j[k]])) {
2611: *mat_j++ = tcol - 1;
2612: *mat_a++ = aa[k];
2613: (*mat_ilen)++;
2615: }
2616: }
2617: }
2618: MatSeqAIJRestoreArrayRead(A,&aa);
2619: /* Free work space */
2620: ISRestoreIndices(iscol,&icol);
2621: PetscFree(smap);
2622: PetscFree(lens);
2623: /* sort */
2624: for (i = 0; i < nrows; i++) {
2625: PetscInt ilen;
2627: mat_i = c->i[i];
2628: mat_j = c->j + mat_i;
2629: mat_a = c->a + mat_i;
2630: ilen = c->ilen[i];
2631: PetscSortIntWithScalarArray(ilen,mat_j,mat_a);
2632: }
2633: }
2634: #if defined(PETSC_HAVE_DEVICE)
2635: MatBindToCPU(C,A->boundtocpu);
2636: #endif
2637: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2638: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2640: ISRestoreIndices(isrow,&irow);
2641: *B = C;
2642: return 0;
2643: }
2645: PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat)
2646: {
2647: Mat B;
2649: if (scall == MAT_INITIAL_MATRIX) {
2650: MatCreate(subComm,&B);
2651: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);
2652: MatSetBlockSizesFromMats(B,mat,mat);
2653: MatSetType(B,MATSEQAIJ);
2654: MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);
2655: *subMat = B;
2656: } else {
2657: MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);
2658: }
2659: return 0;
2660: }
2662: PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2663: {
2664: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2665: Mat outA;
2666: PetscBool row_identity,col_identity;
2670: ISIdentity(row,&row_identity);
2671: ISIdentity(col,&col_identity);
2673: outA = inA;
2674: outA->factortype = MAT_FACTOR_LU;
2675: PetscFree(inA->solvertype);
2676: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
2678: PetscObjectReference((PetscObject)row);
2679: ISDestroy(&a->row);
2681: a->row = row;
2683: PetscObjectReference((PetscObject)col);
2684: ISDestroy(&a->col);
2686: a->col = col;
2688: /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */
2689: ISDestroy(&a->icol);
2690: ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2691: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
2693: if (!a->solve_work) { /* this matrix may have been factored before */
2694: PetscMalloc1(inA->rmap->n+1,&a->solve_work);
2695: PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));
2696: }
2698: MatMarkDiagonal_SeqAIJ(inA);
2699: if (row_identity && col_identity) {
2700: MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);
2701: } else {
2702: MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);
2703: }
2704: return 0;
2705: }
2707: PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha)
2708: {
2709: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2710: PetscScalar *v;
2711: PetscBLASInt one = 1,bnz;
2713: MatSeqAIJGetArray(inA,&v);
2714: PetscBLASIntCast(a->nz,&bnz);
2715: PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one));
2716: PetscLogFlops(a->nz);
2717: MatSeqAIJRestoreArray(inA,&v);
2718: MatSeqAIJInvalidateDiagonal(inA);
2719: return 0;
2720: }
2722: PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj)
2723: {
2724: PetscInt i;
2726: if (!submatj->id) { /* delete data that are linked only to submats[id=0] */
2727: PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);
2729: for (i=0; i<submatj->nrqr; ++i) {
2730: PetscFree(submatj->sbuf2[i]);
2731: }
2732: PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);
2734: if (submatj->rbuf1) {
2735: PetscFree(submatj->rbuf1[0]);
2736: PetscFree(submatj->rbuf1);
2737: }
2739: for (i=0; i<submatj->nrqs; ++i) {
2740: PetscFree(submatj->rbuf3[i]);
2741: }
2742: PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);
2743: PetscFree(submatj->pa);
2744: }
2746: #if defined(PETSC_USE_CTABLE)
2747: PetscTableDestroy((PetscTable*)&submatj->rmap);
2748: if (submatj->cmap_loc) PetscFree(submatj->cmap_loc);
2749: PetscFree(submatj->rmap_loc);
2750: #else
2751: PetscFree(submatj->rmap);
2752: #endif
2754: if (!submatj->allcolumns) {
2755: #if defined(PETSC_USE_CTABLE)
2756: PetscTableDestroy((PetscTable*)&submatj->cmap);
2757: #else
2758: PetscFree(submatj->cmap);
2759: #endif
2760: }
2761: PetscFree(submatj->row2proc);
2763: PetscFree(submatj);
2764: return 0;
2765: }
2767: PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C)
2768: {
2769: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
2770: Mat_SubSppt *submatj = c->submatis1;
2772: (*submatj->destroy)(C);
2773: MatDestroySubMatrix_Private(submatj);
2774: return 0;
2775: }
2777: PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[])
2778: {
2779: PetscInt i;
2780: Mat C;
2781: Mat_SeqAIJ *c;
2782: Mat_SubSppt *submatj;
2784: for (i=0; i<n; i++) {
2785: C = (*mat)[i];
2786: c = (Mat_SeqAIJ*)C->data;
2787: submatj = c->submatis1;
2788: if (submatj) {
2789: if (--((PetscObject)C)->refct <= 0) {
2790: (*submatj->destroy)(C);
2791: MatDestroySubMatrix_Private(submatj);
2792: PetscFree(C->defaultvectype);
2793: PetscLayoutDestroy(&C->rmap);
2794: PetscLayoutDestroy(&C->cmap);
2795: PetscHeaderDestroy(&C);
2796: }
2797: } else {
2798: MatDestroy(&C);
2799: }
2800: }
2802: /* Destroy Dummy submatrices created for reuse */
2803: MatDestroySubMatrices_Dummy(n,mat);
2805: PetscFree(*mat);
2806: return 0;
2807: }
2809: PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2810: {
2811: PetscInt i;
2813: if (scall == MAT_INITIAL_MATRIX) {
2814: PetscCalloc1(n+1,B);
2815: }
2817: for (i=0; i<n; i++) {
2818: MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2819: }
2820: return 0;
2821: }
2823: PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov)
2824: {
2825: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2826: PetscInt row,i,j,k,l,m,n,*nidx,isz,val;
2827: const PetscInt *idx;
2828: PetscInt start,end,*ai,*aj;
2829: PetscBT table;
2831: m = A->rmap->n;
2832: ai = a->i;
2833: aj = a->j;
2837: PetscMalloc1(m+1,&nidx);
2838: PetscBTCreate(m,&table);
2840: for (i=0; i<is_max; i++) {
2841: /* Initialize the two local arrays */
2842: isz = 0;
2843: PetscBTMemzero(m,table);
2845: /* Extract the indices, assume there can be duplicate entries */
2846: ISGetIndices(is[i],&idx);
2847: ISGetLocalSize(is[i],&n);
2849: /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */
2850: for (j=0; j<n; ++j) {
2851: if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j];
2852: }
2853: ISRestoreIndices(is[i],&idx);
2854: ISDestroy(&is[i]);
2856: k = 0;
2857: for (j=0; j<ov; j++) { /* for each overlap */
2858: n = isz;
2859: for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */
2860: row = nidx[k];
2861: start = ai[row];
2862: end = ai[row+1];
2863: for (l = start; l<end; l++) {
2864: val = aj[l];
2865: if (!PetscBTLookupSet(table,val)) nidx[isz++] = val;
2866: }
2867: }
2868: }
2869: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));
2870: }
2871: PetscBTDestroy(&table);
2872: PetscFree(nidx);
2873: return 0;
2874: }
2876: /* -------------------------------------------------------------- */
2877: PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B)
2878: {
2879: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2880: PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n;
2881: const PetscInt *row,*col;
2882: PetscInt *cnew,j,*lens;
2883: IS icolp,irowp;
2884: PetscInt *cwork = NULL;
2885: PetscScalar *vwork = NULL;
2887: ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);
2888: ISGetIndices(irowp,&row);
2889: ISInvertPermutation(colp,PETSC_DECIDE,&icolp);
2890: ISGetIndices(icolp,&col);
2892: /* determine lengths of permuted rows */
2893: PetscMalloc1(m+1,&lens);
2894: for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i];
2895: MatCreate(PetscObjectComm((PetscObject)A),B);
2896: MatSetSizes(*B,m,n,m,n);
2897: MatSetBlockSizesFromMats(*B,A,A);
2898: MatSetType(*B,((PetscObject)A)->type_name);
2899: MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);
2900: PetscFree(lens);
2902: PetscMalloc1(n,&cnew);
2903: for (i=0; i<m; i++) {
2904: MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2905: for (j=0; j<nz; j++) cnew[j] = col[cwork[j]];
2906: MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);
2907: MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2908: }
2909: PetscFree(cnew);
2911: (*B)->assembled = PETSC_FALSE;
2913: #if defined(PETSC_HAVE_DEVICE)
2914: MatBindToCPU(*B,A->boundtocpu);
2915: #endif
2916: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
2917: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
2918: ISRestoreIndices(irowp,&row);
2919: ISRestoreIndices(icolp,&col);
2920: ISDestroy(&irowp);
2921: ISDestroy(&icolp);
2922: if (rowp == colp) {
2923: MatPropagateSymmetryOptions(A,*B);
2924: }
2925: return 0;
2926: }
2928: PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str)
2929: {
2930: /* If the two matrices have the same copy implementation, use fast copy. */
2931: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2932: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2933: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
2934: const PetscScalar *aa;
2936: MatSeqAIJGetArrayRead(A,&aa);
2938: PetscArraycpy(b->a,aa,a->i[A->rmap->n]);
2939: PetscObjectStateIncrease((PetscObject)B);
2940: MatSeqAIJRestoreArrayRead(A,&aa);
2941: } else {
2942: MatCopy_Basic(A,B,str);
2943: }
2944: return 0;
2945: }
2947: PetscErrorCode MatSetUp_SeqAIJ(Mat A)
2948: {
2949: MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL);
2950: return 0;
2951: }
2953: PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[])
2954: {
2955: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2957: *array = a->a;
2958: return 0;
2959: }
2961: PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[])
2962: {
2963: *array = NULL;
2964: return 0;
2965: }
2967: /*
2968: Computes the number of nonzeros per row needed for preallocation when X and Y
2969: have different nonzero structure.
2970: */
2971: PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz)
2972: {
2973: PetscInt i,j,k,nzx,nzy;
2975: /* Set the number of nonzeros in the new matrix */
2976: for (i=0; i<m; i++) {
2977: const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i];
2978: nzx = xi[i+1] - xi[i];
2979: nzy = yi[i+1] - yi[i];
2980: nnz[i] = 0;
2981: for (j=0,k=0; j<nzx; j++) { /* Point in X */
2982: for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */
2983: if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */
2984: nnz[i]++;
2985: }
2986: for (; k<nzy; k++) nnz[i]++;
2987: }
2988: return 0;
2989: }
2991: PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz)
2992: {
2993: PetscInt m = Y->rmap->N;
2994: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data;
2995: Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data;
2997: /* Set the number of nonzeros in the new matrix */
2998: MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);
2999: return 0;
3000: }
3002: PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
3003: {
3004: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data;
3006: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
3007: PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE;
3008: if (e) {
3009: PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e);
3010: if (e) {
3011: PetscArraycmp(x->j,y->j,y->nz,&e);
3012: if (e) str = SAME_NONZERO_PATTERN;
3013: }
3014: }
3016: }
3017: if (str == SAME_NONZERO_PATTERN) {
3018: const PetscScalar *xa;
3019: PetscScalar *ya,alpha = a;
3020: PetscBLASInt one = 1,bnz;
3022: PetscBLASIntCast(x->nz,&bnz);
3023: MatSeqAIJGetArray(Y,&ya);
3024: MatSeqAIJGetArrayRead(X,&xa);
3025: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one));
3026: MatSeqAIJRestoreArrayRead(X,&xa);
3027: MatSeqAIJRestoreArray(Y,&ya);
3028: PetscLogFlops(2.0*bnz);
3029: MatSeqAIJInvalidateDiagonal(Y);
3030: PetscObjectStateIncrease((PetscObject)Y);
3031: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
3032: MatAXPY_Basic(Y,a,X,str);
3033: } else {
3034: Mat B;
3035: PetscInt *nnz;
3036: PetscMalloc1(Y->rmap->N,&nnz);
3037: MatCreate(PetscObjectComm((PetscObject)Y),&B);
3038: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
3039: MatSetLayouts(B,Y->rmap,Y->cmap);
3040: MatSetType(B,((PetscObject)Y)->type_name);
3041: MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);
3042: MatSeqAIJSetPreallocation(B,0,nnz);
3043: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
3044: MatHeaderMerge(Y,&B);
3045: PetscFree(nnz);
3046: }
3047: return 0;
3048: }
3050: PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat)
3051: {
3052: #if defined(PETSC_USE_COMPLEX)
3053: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3054: PetscInt i,nz;
3055: PetscScalar *a;
3057: nz = aij->nz;
3058: MatSeqAIJGetArray(mat,&a);
3059: for (i=0; i<nz; i++) a[i] = PetscConj(a[i]);
3060: MatSeqAIJRestoreArray(mat,&a);
3061: #else
3062: #endif
3063: return 0;
3064: }
3066: PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3067: {
3068: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3069: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3070: PetscReal atmp;
3071: PetscScalar *x;
3072: const MatScalar *aa,*av;
3075: MatSeqAIJGetArrayRead(A,&av);
3076: aa = av;
3077: ai = a->i;
3078: aj = a->j;
3080: VecSet(v,0.0);
3081: VecGetArrayWrite(v,&x);
3082: VecGetLocalSize(v,&n);
3084: for (i=0; i<m; i++) {
3085: ncols = ai[1] - ai[0]; ai++;
3086: for (j=0; j<ncols; j++) {
3087: atmp = PetscAbsScalar(*aa);
3088: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
3089: aa++; aj++;
3090: }
3091: }
3092: VecRestoreArrayWrite(v,&x);
3093: MatSeqAIJRestoreArrayRead(A,&av);
3094: return 0;
3095: }
3097: PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3098: {
3099: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3100: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3101: PetscScalar *x;
3102: const MatScalar *aa,*av;
3105: MatSeqAIJGetArrayRead(A,&av);
3106: aa = av;
3107: ai = a->i;
3108: aj = a->j;
3110: VecSet(v,0.0);
3111: VecGetArrayWrite(v,&x);
3112: VecGetLocalSize(v,&n);
3114: for (i=0; i<m; i++) {
3115: ncols = ai[1] - ai[0]; ai++;
3116: if (ncols == A->cmap->n) { /* row is dense */
3117: x[i] = *aa; if (idx) idx[i] = 0;
3118: } else { /* row is sparse so already KNOW maximum is 0.0 or higher */
3119: x[i] = 0.0;
3120: if (idx) {
3121: for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */
3122: if (aj[j] > j) {
3123: idx[i] = j;
3124: break;
3125: }
3126: }
3127: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3128: if (j==ncols && j < A->cmap->n) idx[i] = j;
3129: }
3130: }
3131: for (j=0; j<ncols; j++) {
3132: if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3133: aa++; aj++;
3134: }
3135: }
3136: VecRestoreArrayWrite(v,&x);
3137: MatSeqAIJRestoreArrayRead(A,&av);
3138: return 0;
3139: }
3141: PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3142: {
3143: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3144: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3145: PetscScalar *x;
3146: const MatScalar *aa,*av;
3148: MatSeqAIJGetArrayRead(A,&av);
3149: aa = av;
3150: ai = a->i;
3151: aj = a->j;
3153: VecSet(v,0.0);
3154: VecGetArrayWrite(v,&x);
3155: VecGetLocalSize(v,&n);
3157: for (i=0; i<m; i++) {
3158: ncols = ai[1] - ai[0]; ai++;
3159: if (ncols == A->cmap->n) { /* row is dense */
3160: x[i] = *aa; if (idx) idx[i] = 0;
3161: } else { /* row is sparse so already KNOW minimum is 0.0 or higher */
3162: x[i] = 0.0;
3163: if (idx) { /* find first implicit 0.0 in the row */
3164: for (j=0; j<ncols; j++) {
3165: if (aj[j] > j) {
3166: idx[i] = j;
3167: break;
3168: }
3169: }
3170: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3171: if (j==ncols && j < A->cmap->n) idx[i] = j;
3172: }
3173: }
3174: for (j=0; j<ncols; j++) {
3175: if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3176: aa++; aj++;
3177: }
3178: }
3179: VecRestoreArrayWrite(v,&x);
3180: MatSeqAIJRestoreArrayRead(A,&av);
3181: return 0;
3182: }
3184: PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3185: {
3186: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3187: PetscInt i,j,m = A->rmap->n,ncols,n;
3188: const PetscInt *ai,*aj;
3189: PetscScalar *x;
3190: const MatScalar *aa,*av;
3193: MatSeqAIJGetArrayRead(A,&av);
3194: aa = av;
3195: ai = a->i;
3196: aj = a->j;
3198: VecSet(v,0.0);
3199: VecGetArrayWrite(v,&x);
3200: VecGetLocalSize(v,&n);
3202: for (i=0; i<m; i++) {
3203: ncols = ai[1] - ai[0]; ai++;
3204: if (ncols == A->cmap->n) { /* row is dense */
3205: x[i] = *aa; if (idx) idx[i] = 0;
3206: } else { /* row is sparse so already KNOW minimum is 0.0 or lower */
3207: x[i] = 0.0;
3208: if (idx) { /* find first implicit 0.0 in the row */
3209: for (j=0; j<ncols; j++) {
3210: if (aj[j] > j) {
3211: idx[i] = j;
3212: break;
3213: }
3214: }
3215: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3216: if (j==ncols && j < A->cmap->n) idx[i] = j;
3217: }
3218: }
3219: for (j=0; j<ncols; j++) {
3220: if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3221: aa++; aj++;
3222: }
3223: }
3224: VecRestoreArrayWrite(v,&x);
3225: MatSeqAIJRestoreArrayRead(A,&av);
3226: return 0;
3227: }
3229: PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values)
3230: {
3231: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
3232: PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j;
3233: MatScalar *diag,work[25],*v_work;
3234: const PetscReal shift = 0.0;
3235: PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
3237: allowzeropivot = PetscNot(A->erroriffailure);
3238: if (a->ibdiagvalid) {
3239: if (values) *values = a->ibdiag;
3240: return 0;
3241: }
3242: MatMarkDiagonal_SeqAIJ(A);
3243: if (!a->ibdiag) {
3244: PetscMalloc1(bs2*mbs,&a->ibdiag);
3245: PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));
3246: }
3247: diag = a->ibdiag;
3248: if (values) *values = a->ibdiag;
3249: /* factor and invert each block */
3250: switch (bs) {
3251: case 1:
3252: for (i=0; i<mbs; i++) {
3253: MatGetValues(A,1,&i,1,&i,diag+i);
3254: if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) {
3255: if (allowzeropivot) {
3256: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3257: A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]);
3258: A->factorerror_zeropivot_row = i;
3259: PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);
3260: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);
3261: }
3262: diag[i] = (PetscScalar)1.0 / (diag[i] + shift);
3263: }
3264: break;
3265: case 2:
3266: for (i=0; i<mbs; i++) {
3267: ij[0] = 2*i; ij[1] = 2*i + 1;
3268: MatGetValues(A,2,ij,2,ij,diag);
3269: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
3270: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3271: PetscKernel_A_gets_transpose_A_2(diag);
3272: diag += 4;
3273: }
3274: break;
3275: case 3:
3276: for (i=0; i<mbs; i++) {
3277: ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2;
3278: MatGetValues(A,3,ij,3,ij,diag);
3279: PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
3280: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3281: PetscKernel_A_gets_transpose_A_3(diag);
3282: diag += 9;
3283: }
3284: break;
3285: case 4:
3286: for (i=0; i<mbs; i++) {
3287: ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3;
3288: MatGetValues(A,4,ij,4,ij,diag);
3289: PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
3290: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3291: PetscKernel_A_gets_transpose_A_4(diag);
3292: diag += 16;
3293: }
3294: break;
3295: case 5:
3296: for (i=0; i<mbs; i++) {
3297: ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4;
3298: MatGetValues(A,5,ij,5,ij,diag);
3299: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
3300: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3301: PetscKernel_A_gets_transpose_A_5(diag);
3302: diag += 25;
3303: }
3304: break;
3305: case 6:
3306: for (i=0; i<mbs; i++) {
3307: ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5;
3308: MatGetValues(A,6,ij,6,ij,diag);
3309: PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
3310: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3311: PetscKernel_A_gets_transpose_A_6(diag);
3312: diag += 36;
3313: }
3314: break;
3315: case 7:
3316: for (i=0; i<mbs; i++) {
3317: ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6;
3318: MatGetValues(A,7,ij,7,ij,diag);
3319: PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
3320: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3321: PetscKernel_A_gets_transpose_A_7(diag);
3322: diag += 49;
3323: }
3324: break;
3325: default:
3326: PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);
3327: for (i=0; i<mbs; i++) {
3328: for (j=0; j<bs; j++) {
3329: IJ[j] = bs*i + j;
3330: }
3331: MatGetValues(A,bs,IJ,bs,IJ,diag);
3332: PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
3333: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3334: PetscKernel_A_gets_transpose_A_N(diag,bs);
3335: diag += bs2;
3336: }
3337: PetscFree3(v_work,v_pivots,IJ);
3338: }
3339: a->ibdiagvalid = PETSC_TRUE;
3340: return 0;
3341: }
3343: static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx)
3344: {
3345: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3346: PetscScalar a,*aa;
3347: PetscInt m,n,i,j,col;
3349: if (!x->assembled) {
3350: MatGetSize(x,&m,&n);
3351: for (i=0; i<m; i++) {
3352: for (j=0; j<aij->imax[i]; j++) {
3353: PetscRandomGetValue(rctx,&a);
3354: col = (PetscInt)(n*PetscRealPart(a));
3355: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3356: }
3357: }
3358: } else {
3359: MatSeqAIJGetArrayWrite(x,&aa);
3360: for (i=0; i<aij->nz; i++) PetscRandomGetValue(rctx,aa+i);
3361: MatSeqAIJRestoreArrayWrite(x,&aa);
3362: }
3363: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3364: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3365: return 0;
3366: }
3368: /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */
3369: PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx)
3370: {
3371: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3372: PetscScalar a;
3373: PetscInt m,n,i,j,col,nskip;
3375: nskip = high - low;
3376: MatGetSize(x,&m,&n);
3377: n -= nskip; /* shrink number of columns where nonzeros can be set */
3378: for (i=0; i<m; i++) {
3379: for (j=0; j<aij->imax[i]; j++) {
3380: PetscRandomGetValue(rctx,&a);
3381: col = (PetscInt)(n*PetscRealPart(a));
3382: if (col >= low) col += nskip; /* shift col rightward to skip the hole */
3383: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3384: }
3385: }
3386: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3387: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3388: return 0;
3389: }
3391: /* -------------------------------------------------------------------*/
3392: static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ,
3393: MatGetRow_SeqAIJ,
3394: MatRestoreRow_SeqAIJ,
3395: MatMult_SeqAIJ,
3396: /* 4*/ MatMultAdd_SeqAIJ,
3397: MatMultTranspose_SeqAIJ,
3398: MatMultTransposeAdd_SeqAIJ,
3399: NULL,
3400: NULL,
3401: NULL,
3402: /* 10*/ NULL,
3403: MatLUFactor_SeqAIJ,
3404: NULL,
3405: MatSOR_SeqAIJ,
3406: MatTranspose_SeqAIJ,
3407: /*1 5*/ MatGetInfo_SeqAIJ,
3408: MatEqual_SeqAIJ,
3409: MatGetDiagonal_SeqAIJ,
3410: MatDiagonalScale_SeqAIJ,
3411: MatNorm_SeqAIJ,
3412: /* 20*/ NULL,
3413: MatAssemblyEnd_SeqAIJ,
3414: MatSetOption_SeqAIJ,
3415: MatZeroEntries_SeqAIJ,
3416: /* 24*/ MatZeroRows_SeqAIJ,
3417: NULL,
3418: NULL,
3419: NULL,
3420: NULL,
3421: /* 29*/ MatSetUp_SeqAIJ,
3422: NULL,
3423: NULL,
3424: NULL,
3425: NULL,
3426: /* 34*/ MatDuplicate_SeqAIJ,
3427: NULL,
3428: NULL,
3429: MatILUFactor_SeqAIJ,
3430: NULL,
3431: /* 39*/ MatAXPY_SeqAIJ,
3432: MatCreateSubMatrices_SeqAIJ,
3433: MatIncreaseOverlap_SeqAIJ,
3434: MatGetValues_SeqAIJ,
3435: MatCopy_SeqAIJ,
3436: /* 44*/ MatGetRowMax_SeqAIJ,
3437: MatScale_SeqAIJ,
3438: MatShift_SeqAIJ,
3439: MatDiagonalSet_SeqAIJ,
3440: MatZeroRowsColumns_SeqAIJ,
3441: /* 49*/ MatSetRandom_SeqAIJ,
3442: MatGetRowIJ_SeqAIJ,
3443: MatRestoreRowIJ_SeqAIJ,
3444: MatGetColumnIJ_SeqAIJ,
3445: MatRestoreColumnIJ_SeqAIJ,
3446: /* 54*/ MatFDColoringCreate_SeqXAIJ,
3447: NULL,
3448: NULL,
3449: MatPermute_SeqAIJ,
3450: NULL,
3451: /* 59*/ NULL,
3452: MatDestroy_SeqAIJ,
3453: MatView_SeqAIJ,
3454: NULL,
3455: NULL,
3456: /* 64*/ NULL,
3457: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ,
3458: NULL,
3459: NULL,
3460: NULL,
3461: /* 69*/ MatGetRowMaxAbs_SeqAIJ,
3462: MatGetRowMinAbs_SeqAIJ,
3463: NULL,
3464: NULL,
3465: NULL,
3466: /* 74*/ NULL,
3467: MatFDColoringApply_AIJ,
3468: NULL,
3469: NULL,
3470: NULL,
3471: /* 79*/ MatFindZeroDiagonals_SeqAIJ,
3472: NULL,
3473: NULL,
3474: NULL,
3475: MatLoad_SeqAIJ,
3476: /* 84*/ MatIsSymmetric_SeqAIJ,
3477: MatIsHermitian_SeqAIJ,
3478: NULL,
3479: NULL,
3480: NULL,
3481: /* 89*/ NULL,
3482: NULL,
3483: MatMatMultNumeric_SeqAIJ_SeqAIJ,
3484: NULL,
3485: NULL,
3486: /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy,
3487: NULL,
3488: NULL,
3489: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ,
3490: NULL,
3491: /* 99*/ MatProductSetFromOptions_SeqAIJ,
3492: NULL,
3493: NULL,
3494: MatConjugate_SeqAIJ,
3495: NULL,
3496: /*104*/ MatSetValuesRow_SeqAIJ,
3497: MatRealPart_SeqAIJ,
3498: MatImaginaryPart_SeqAIJ,
3499: NULL,
3500: NULL,
3501: /*109*/ MatMatSolve_SeqAIJ,
3502: NULL,
3503: MatGetRowMin_SeqAIJ,
3504: NULL,
3505: MatMissingDiagonal_SeqAIJ,
3506: /*114*/ NULL,
3507: NULL,
3508: NULL,
3509: NULL,
3510: NULL,
3511: /*119*/ NULL,
3512: NULL,
3513: NULL,
3514: NULL,
3515: MatGetMultiProcBlock_SeqAIJ,
3516: /*124*/ MatFindNonzeroRows_SeqAIJ,
3517: MatGetColumnReductions_SeqAIJ,
3518: MatInvertBlockDiagonal_SeqAIJ,
3519: MatInvertVariableBlockDiagonal_SeqAIJ,
3520: NULL,
3521: /*129*/ NULL,
3522: NULL,
3523: NULL,
3524: MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ,
3525: MatTransposeColoringCreate_SeqAIJ,
3526: /*134*/ MatTransColoringApplySpToDen_SeqAIJ,
3527: MatTransColoringApplyDenToSp_SeqAIJ,
3528: NULL,
3529: NULL,
3530: MatRARtNumeric_SeqAIJ_SeqAIJ,
3531: /*139*/NULL,
3532: NULL,
3533: NULL,
3534: MatFDColoringSetUp_SeqXAIJ,
3535: MatFindOffBlockDiagonalEntries_SeqAIJ,
3536: MatCreateMPIMatConcatenateSeqMat_SeqAIJ,
3537: /*145*/MatDestroySubMatrices_SeqAIJ,
3538: NULL,
3539: NULL
3540: };
3542: PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices)
3543: {
3544: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3545: PetscInt i,nz,n;
3547: nz = aij->maxnz;
3548: n = mat->rmap->n;
3549: for (i=0; i<nz; i++) {
3550: aij->j[i] = indices[i];
3551: }
3552: aij->nz = nz;
3553: for (i=0; i<n; i++) {
3554: aij->ilen[i] = aij->imax[i];
3555: }
3556: return 0;
3557: }
3559: /*
3560: * Given a sparse matrix with global column indices, compact it by using a local column space.
3561: * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable()
3562: */
3563: PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping)
3564: {
3565: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3566: PetscTable gid1_lid1;
3567: PetscTablePosition tpos;
3568: PetscInt gid,lid,i,ec,nz = aij->nz;
3569: PetscInt *garray,*jj = aij->j;
3573: /* use a table */
3574: PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);
3575: ec = 0;
3576: for (i=0; i<nz; i++) {
3577: PetscInt data,gid1 = jj[i] + 1;
3578: PetscTableFind(gid1_lid1,gid1,&data);
3579: if (!data) {
3580: /* one based table */
3581: PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);
3582: }
3583: }
3584: /* form array of columns we need */
3585: PetscMalloc1(ec,&garray);
3586: PetscTableGetHeadPosition(gid1_lid1,&tpos);
3587: while (tpos) {
3588: PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);
3589: gid--;
3590: lid--;
3591: garray[lid] = gid;
3592: }
3593: PetscSortInt(ec,garray); /* sort, and rebuild */
3594: PetscTableRemoveAll(gid1_lid1);
3595: for (i=0; i<ec; i++) {
3596: PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);
3597: }
3598: /* compact out the extra columns in B */
3599: for (i=0; i<nz; i++) {
3600: PetscInt gid1 = jj[i] + 1;
3601: PetscTableFind(gid1_lid1,gid1,&lid);
3602: lid--;
3603: jj[i] = lid;
3604: }
3605: PetscLayoutDestroy(&mat->cmap);
3606: PetscTableDestroy(&gid1_lid1);
3607: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap);
3608: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);
3609: ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);
3610: return 0;
3611: }
3613: /*@
3614: MatSeqAIJSetColumnIndices - Set the column indices for all the rows
3615: in the matrix.
3617: Input Parameters:
3618: + mat - the SeqAIJ matrix
3619: - indices - the column indices
3621: Level: advanced
3623: Notes:
3624: This can be called if you have precomputed the nonzero structure of the
3625: matrix and want to provide it to the matrix object to improve the performance
3626: of the MatSetValues() operation.
3628: You MUST have set the correct numbers of nonzeros per row in the call to
3629: MatCreateSeqAIJ(), and the columns indices MUST be sorted.
3631: MUST be called before any calls to MatSetValues();
3633: The indices should start with zero, not one.
3635: @*/
3636: PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices)
3637: {
3640: PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
3641: return 0;
3642: }
3644: /* ----------------------------------------------------------------------------------------*/
3646: PetscErrorCode MatStoreValues_SeqAIJ(Mat mat)
3647: {
3648: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3649: size_t nz = aij->i[mat->rmap->n];
3653: /* allocate space for values if not already there */
3654: if (!aij->saved_values) {
3655: PetscMalloc1(nz+1,&aij->saved_values);
3656: PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
3657: }
3659: /* copy values over */
3660: PetscArraycpy(aij->saved_values,aij->a,nz);
3661: return 0;
3662: }
3664: /*@
3665: MatStoreValues - Stashes a copy of the matrix values; this allows, for
3666: example, reuse of the linear part of a Jacobian, while recomputing the
3667: nonlinear portion.
3669: Collect on Mat
3671: Input Parameters:
3672: . mat - the matrix (currently only AIJ matrices support this option)
3674: Level: advanced
3676: Common Usage, with SNESSolve():
3677: $ Create Jacobian matrix
3678: $ Set linear terms into matrix
3679: $ Apply boundary conditions to matrix, at this time matrix must have
3680: $ final nonzero structure (i.e. setting the nonlinear terms and applying
3681: $ boundary conditions again will not change the nonzero structure
3682: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3683: $ MatStoreValues(mat);
3684: $ Call SNESSetJacobian() with matrix
3685: $ In your Jacobian routine
3686: $ MatRetrieveValues(mat);
3687: $ Set nonlinear terms in matrix
3689: Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself:
3690: $ // build linear portion of Jacobian
3691: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3692: $ MatStoreValues(mat);
3693: $ loop over nonlinear iterations
3694: $ MatRetrieveValues(mat);
3695: $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian
3696: $ // call MatAssemblyBegin/End() on matrix
3697: $ Solve linear system with Jacobian
3698: $ endloop
3700: Notes:
3701: Matrix must already be assemblied before calling this routine
3702: Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before
3703: calling this routine.
3705: When this is called multiple times it overwrites the previous set of stored values
3706: and does not allocated additional space.
3708: .seealso: MatRetrieveValues()
3710: @*/
3711: PetscErrorCode MatStoreValues(Mat mat)
3712: {
3716: PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));
3717: return 0;
3718: }
3720: PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat)
3721: {
3722: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3723: PetscInt nz = aij->i[mat->rmap->n];
3727: /* copy values over */
3728: PetscArraycpy(aij->a,aij->saved_values,nz);
3729: return 0;
3730: }
3732: /*@
3733: MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for
3734: example, reuse of the linear part of a Jacobian, while recomputing the
3735: nonlinear portion.
3737: Collect on Mat
3739: Input Parameters:
3740: . mat - the matrix (currently only AIJ matrices support this option)
3742: Level: advanced
3744: .seealso: MatStoreValues()
3746: @*/
3747: PetscErrorCode MatRetrieveValues(Mat mat)
3748: {
3752: PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));
3753: return 0;
3754: }
3756: /* --------------------------------------------------------------------------------*/
3757: /*@C
3758: MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format
3759: (the default parallel PETSc format). For good matrix assembly performance
3760: the user should preallocate the matrix storage by setting the parameter nz
3761: (or the array nnz). By setting these parameters accurately, performance
3762: during matrix assembly can be increased by more than a factor of 50.
3764: Collective
3766: Input Parameters:
3767: + comm - MPI communicator, set to PETSC_COMM_SELF
3768: . m - number of rows
3769: . n - number of columns
3770: . nz - number of nonzeros per row (same for all rows)
3771: - nnz - array containing the number of nonzeros in the various rows
3772: (possibly different for each row) or NULL
3774: Output Parameter:
3775: . A - the matrix
3777: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3778: MatXXXXSetPreallocation() paradigm instead of this routine directly.
3779: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
3781: Notes:
3782: If nnz is given then nz is ignored
3784: The AIJ format (also called the Yale sparse matrix format or
3785: compressed row storage), is fully compatible with standard Fortran 77
3786: storage. That is, the stored row and column indices can begin at
3787: either one (as in Fortran) or zero. See the users' manual for details.
3789: Specify the preallocated storage with either nz or nnz (not both).
3790: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3791: allocation. For large problems you MUST preallocate memory or you
3792: will get TERRIBLE performance, see the users' manual chapter on matrices.
3794: By default, this format uses inodes (identical nodes) when possible, to
3795: improve numerical efficiency of matrix-vector products and solves. We
3796: search for consecutive rows with the same nonzero structure, thereby
3797: reusing matrix information to achieve increased efficiency.
3799: Options Database Keys:
3800: + -mat_no_inode - Do not use inodes
3801: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3803: Level: intermediate
3805: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays()
3807: @*/
3808: PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3809: {
3810: MatCreate(comm,A);
3811: MatSetSizes(*A,m,n,m,n);
3812: MatSetType(*A,MATSEQAIJ);
3813: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
3814: return 0;
3815: }
3817: /*@C
3818: MatSeqAIJSetPreallocation - For good matrix assembly performance
3819: the user should preallocate the matrix storage by setting the parameter nz
3820: (or the array nnz). By setting these parameters accurately, performance
3821: during matrix assembly can be increased by more than a factor of 50.
3823: Collective
3825: Input Parameters:
3826: + B - The matrix
3827: . nz - number of nonzeros per row (same for all rows)
3828: - nnz - array containing the number of nonzeros in the various rows
3829: (possibly different for each row) or NULL
3831: Notes:
3832: If nnz is given then nz is ignored
3834: The AIJ format (also called the Yale sparse matrix format or
3835: compressed row storage), is fully compatible with standard Fortran 77
3836: storage. That is, the stored row and column indices can begin at
3837: either one (as in Fortran) or zero. See the users' manual for details.
3839: Specify the preallocated storage with either nz or nnz (not both).
3840: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3841: allocation. For large problems you MUST preallocate memory or you
3842: will get TERRIBLE performance, see the users' manual chapter on matrices.
3844: You can call MatGetInfo() to get information on how effective the preallocation was;
3845: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3846: You can also run with the option -info and look for messages with the string
3847: malloc in them to see if additional memory allocation was needed.
3849: Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix
3850: entries or columns indices
3852: By default, this format uses inodes (identical nodes) when possible, to
3853: improve numerical efficiency of matrix-vector products and solves. We
3854: search for consecutive rows with the same nonzero structure, thereby
3855: reusing matrix information to achieve increased efficiency.
3857: Options Database Keys:
3858: + -mat_no_inode - Do not use inodes
3859: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3861: Level: intermediate
3863: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo(),
3864: MatSeqAIJSetTotalPreallocation()
3866: @*/
3867: PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[])
3868: {
3871: PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));
3872: return 0;
3873: }
3875: PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz)
3876: {
3877: Mat_SeqAIJ *b;
3878: PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;
3879: PetscInt i;
3881: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3882: if (nz == MAT_SKIP_ALLOCATION) {
3883: skipallocation = PETSC_TRUE;
3884: nz = 0;
3885: }
3886: PetscLayoutSetUp(B->rmap);
3887: PetscLayoutSetUp(B->cmap);
3889: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3891: if (PetscUnlikelyDebug(nnz)) {
3892: for (i=0; i<B->rmap->n; i++) {
3895: }
3896: }
3898: B->preallocated = PETSC_TRUE;
3900: b = (Mat_SeqAIJ*)B->data;
3902: if (!skipallocation) {
3903: if (!b->imax) {
3904: PetscMalloc1(B->rmap->n,&b->imax);
3905: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3906: }
3907: if (!b->ilen) {
3908: /* b->ilen will count nonzeros in each row so far. */
3909: PetscCalloc1(B->rmap->n,&b->ilen);
3910: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3911: } else {
3912: PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt));
3913: }
3914: if (!b->ipre) {
3915: PetscMalloc1(B->rmap->n,&b->ipre);
3916: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3917: }
3918: if (!nnz) {
3919: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10;
3920: else if (nz < 0) nz = 1;
3921: nz = PetscMin(nz,B->cmap->n);
3922: for (i=0; i<B->rmap->n; i++) b->imax[i] = nz;
3923: nz = nz*B->rmap->n;
3924: } else {
3925: PetscInt64 nz64 = 0;
3926: for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
3927: PetscIntCast(nz64,&nz);
3928: }
3930: /* allocate the matrix space */
3931: /* FIXME: should B's old memory be unlogged? */
3932: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
3933: if (B->structure_only) {
3934: PetscMalloc1(nz,&b->j);
3935: PetscMalloc1(B->rmap->n+1,&b->i);
3936: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));
3937: } else {
3938: PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);
3939: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
3940: }
3941: b->i[0] = 0;
3942: for (i=1; i<B->rmap->n+1; i++) {
3943: b->i[i] = b->i[i-1] + b->imax[i-1];
3944: }
3945: if (B->structure_only) {
3946: b->singlemalloc = PETSC_FALSE;
3947: b->free_a = PETSC_FALSE;
3948: } else {
3949: b->singlemalloc = PETSC_TRUE;
3950: b->free_a = PETSC_TRUE;
3951: }
3952: b->free_ij = PETSC_TRUE;
3953: } else {
3954: b->free_a = PETSC_FALSE;
3955: b->free_ij = PETSC_FALSE;
3956: }
3958: if (b->ipre && nnz != b->ipre && b->imax) {
3959: /* reserve user-requested sparsity */
3960: PetscArraycpy(b->ipre,b->imax,B->rmap->n);
3961: }
3963: b->nz = 0;
3964: b->maxnz = nz;
3965: B->info.nz_unneeded = (double)b->maxnz;
3966: if (realalloc) {
3967: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
3968: }
3969: B->was_assembled = PETSC_FALSE;
3970: B->assembled = PETSC_FALSE;
3971: return 0;
3972: }
3974: PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A)
3975: {
3976: Mat_SeqAIJ *a;
3977: PetscInt i;
3981: /* Check local size. If zero, then return */
3982: if (!A->rmap->n) return 0;
3984: a = (Mat_SeqAIJ*)A->data;
3985: /* if no saved info, we error out */
3990: PetscArraycpy(a->imax,a->ipre,A->rmap->n);
3991: PetscArrayzero(a->ilen,A->rmap->n);
3992: a->i[0] = 0;
3993: for (i=1; i<A->rmap->n+1; i++) {
3994: a->i[i] = a->i[i-1] + a->imax[i-1];
3995: }
3996: A->preallocated = PETSC_TRUE;
3997: a->nz = 0;
3998: a->maxnz = a->i[A->rmap->n];
3999: A->info.nz_unneeded = (double)a->maxnz;
4000: A->was_assembled = PETSC_FALSE;
4001: A->assembled = PETSC_FALSE;
4002: return 0;
4003: }
4005: /*@
4006: MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format.
4008: Input Parameters:
4009: + B - the matrix
4010: . i - the indices into j for the start of each row (starts with zero)
4011: . j - the column indices for each row (starts with zero) these must be sorted for each row
4012: - v - optional values in the matrix
4014: Level: developer
4016: Notes:
4017: The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays()
4019: This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero
4020: structure will be the union of all the previous nonzero structures.
4022: Developer Notes:
4023: An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and
4024: then just copies the v values directly with PetscMemcpy().
4026: This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them.
4028: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation()
4029: @*/
4030: PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[])
4031: {
4034: PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));
4035: return 0;
4036: }
4038: PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[])
4039: {
4040: PetscInt i;
4041: PetscInt m,n;
4042: PetscInt nz;
4043: PetscInt *nnz;
4047: PetscLayoutSetUp(B->rmap);
4048: PetscLayoutSetUp(B->cmap);
4050: MatGetSize(B, &m, &n);
4051: PetscMalloc1(m+1, &nnz);
4052: for (i = 0; i < m; i++) {
4053: nz = Ii[i+1]- Ii[i];
4055: nnz[i] = nz;
4056: }
4057: MatSeqAIJSetPreallocation(B, 0, nnz);
4058: PetscFree(nnz);
4060: for (i = 0; i < m; i++) {
4061: MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES);
4062: }
4064: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
4065: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
4067: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
4068: return 0;
4069: }
4071: /*@
4072: MatSeqAIJKron - Computes C, the Kronecker product of A and B.
4074: Input Parameters:
4075: + A - left-hand side matrix
4076: . B - right-hand side matrix
4077: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4079: Output Parameter:
4080: . C - Kronecker product of A and B
4082: Level: intermediate
4084: Notes:
4085: MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron().
4087: .seealso: MatCreateSeqAIJ(), MATSEQAIJ, MATKAIJ, MatReuse
4088: @*/
4089: PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C)
4090: {
4096: if (reuse == MAT_REUSE_MATRIX) {
4099: }
4100: PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C));
4101: return 0;
4102: }
4104: PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C)
4105: {
4106: Mat newmat;
4107: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4108: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
4109: PetscScalar *v;
4110: const PetscScalar *aa,*ba;
4111: PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n;
4112: PetscBool flg;
4118: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg);
4121: if (reuse == MAT_INITIAL_MATRIX) {
4122: PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j);
4123: MatCreate(PETSC_COMM_SELF,&newmat);
4124: MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn);
4125: MatSetType(newmat,MATAIJ);
4126: i[0] = 0;
4127: for (m = 0; m < am; ++m) {
4128: for (p = 0; p < bm; ++p) {
4129: i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]);
4130: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4131: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4132: j[nnz++] = a->j[n]*bn + b->j[q];
4133: }
4134: }
4135: }
4136: }
4137: MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL);
4138: *C = newmat;
4139: PetscFree2(i,j);
4140: nnz = 0;
4141: }
4142: MatSeqAIJGetArray(*C,&v);
4143: MatSeqAIJGetArrayRead(A,&aa);
4144: MatSeqAIJGetArrayRead(B,&ba);
4145: for (m = 0; m < am; ++m) {
4146: for (p = 0; p < bm; ++p) {
4147: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4148: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4149: v[nnz++] = aa[n] * ba[q];
4150: }
4151: }
4152: }
4153: }
4154: MatSeqAIJRestoreArray(*C,&v);
4155: MatSeqAIJRestoreArrayRead(A,&aa);
4156: MatSeqAIJRestoreArrayRead(B,&ba);
4157: return 0;
4158: }
4160: #include <../src/mat/impls/dense/seq/dense.h>
4161: #include <petsc/private/kernels/petscaxpy.h>
4163: /*
4164: Computes (B'*A')' since computing B*A directly is untenable
4166: n p p
4167: [ ] [ ] [ ]
4168: m [ A ] * n [ B ] = m [ C ]
4169: [ ] [ ] [ ]
4171: */
4172: PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C)
4173: {
4174: Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data;
4175: Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data;
4176: Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data;
4177: PetscInt i,j,n,m,q,p;
4178: const PetscInt *ii,*idx;
4179: const PetscScalar *b,*a,*a_q;
4180: PetscScalar *c,*c_q;
4181: PetscInt clda = sub_c->lda;
4182: PetscInt alda = sub_a->lda;
4184: m = A->rmap->n;
4185: n = A->cmap->n;
4186: p = B->cmap->n;
4187: a = sub_a->v;
4188: b = sub_b->a;
4189: c = sub_c->v;
4190: if (clda == m) {
4191: PetscArrayzero(c,m*p);
4192: } else {
4193: for (j=0;j<p;j++)
4194: for (i=0;i<m;i++)
4195: c[j*clda + i] = 0.0;
4196: }
4197: ii = sub_b->i;
4198: idx = sub_b->j;
4199: for (i=0; i<n; i++) {
4200: q = ii[i+1] - ii[i];
4201: while (q-->0) {
4202: c_q = c + clda*(*idx);
4203: a_q = a + alda*i;
4204: PetscKernelAXPY(c_q,*b,a_q,m);
4205: idx++;
4206: b++;
4207: }
4208: }
4209: return 0;
4210: }
4212: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
4213: {
4214: PetscInt m=A->rmap->n,n=B->cmap->n;
4215: PetscBool cisdense;
4218: MatSetSizes(C,m,n,m,n);
4219: MatSetBlockSizesFromMats(C,A,B);
4220: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
4221: if (!cisdense) {
4222: MatSetType(C,MATDENSE);
4223: }
4224: MatSetUp(C);
4226: C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ;
4227: return 0;
4228: }
4230: /* ----------------------------------------------------------------*/
4231: /*MC
4232: MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices,
4233: based on compressed sparse row format.
4235: Options Database Keys:
4236: . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions()
4238: Level: beginner
4240: Notes:
4241: MatSetValues() may be called for this matrix type with a NULL argument for the numerical values,
4242: in this case the values associated with the rows and columns one passes in are set to zero
4243: in the matrix
4245: MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no
4246: space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored
4248: Developer Notes:
4249: It would be nice if all matrix formats supported passing NULL in for the numerical values
4251: .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType, MATSELL, MATSEQSELL, MATMPISELL
4252: M*/
4254: /*MC
4255: MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices.
4257: This matrix type is identical to MATSEQAIJ when constructed with a single process communicator,
4258: and MATMPIAIJ otherwise. As a result, for single process communicators,
4259: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
4260: for communicators controlling multiple processes. It is recommended that you call both of
4261: the above preallocation routines for simplicity.
4263: Options Database Keys:
4264: . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions()
4266: Developer Notes:
4267: Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when
4268: enough exist.
4270: Level: beginner
4272: .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ, MATMPIAIJ, MATSELL, MATSEQSELL, MATMPISELL
4273: M*/
4275: /*MC
4276: MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices.
4278: This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator,
4279: and MATMPIAIJCRL otherwise. As a result, for single process communicators,
4280: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
4281: for communicators controlling multiple processes. It is recommended that you call both of
4282: the above preallocation routines for simplicity.
4284: Options Database Keys:
4285: . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions()
4287: Level: beginner
4289: .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL
4290: M*/
4292: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*);
4293: #if defined(PETSC_HAVE_ELEMENTAL)
4294: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
4295: #endif
4296: #if defined(PETSC_HAVE_SCALAPACK)
4297: PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
4298: #endif
4299: #if defined(PETSC_HAVE_HYPRE)
4300: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*);
4301: #endif
4303: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*);
4304: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*);
4305: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat);
4307: /*@C
4308: MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored
4310: Not Collective
4312: Input Parameter:
4313: . mat - a MATSEQAIJ matrix
4315: Output Parameter:
4316: . array - pointer to the data
4318: Level: intermediate
4320: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4321: @*/
4322: PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array)
4323: {
4324: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4326: if (aij->ops->getarray) {
4327: (*aij->ops->getarray)(A,array);
4328: } else {
4329: *array = aij->a;
4330: }
4331: return 0;
4332: }
4334: /*@C
4335: MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray()
4337: Not Collective
4339: Input Parameters:
4340: + mat - a MATSEQAIJ matrix
4341: - array - pointer to the data
4343: Level: intermediate
4345: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90()
4346: @*/
4347: PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array)
4348: {
4349: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4351: if (aij->ops->restorearray) {
4352: (*aij->ops->restorearray)(A,array);
4353: } else {
4354: *array = NULL;
4355: }
4356: MatSeqAIJInvalidateDiagonal(A);
4357: PetscObjectStateIncrease((PetscObject)A);
4358: return 0;
4359: }
4361: /*@C
4362: MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored
4364: Not Collective
4366: Input Parameter:
4367: . mat - a MATSEQAIJ matrix
4369: Output Parameter:
4370: . array - pointer to the data
4372: Level: intermediate
4374: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead()
4375: @*/
4376: PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array)
4377: {
4378: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4380: if (aij->ops->getarrayread) {
4381: (*aij->ops->getarrayread)(A,array);
4382: } else {
4383: *array = aij->a;
4384: }
4385: return 0;
4386: }
4388: /*@C
4389: MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead
4391: Not Collective
4393: Input Parameter:
4394: . mat - a MATSEQAIJ matrix
4396: Output Parameter:
4397: . array - pointer to the data
4399: Level: intermediate
4401: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4402: @*/
4403: PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array)
4404: {
4405: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4407: if (aij->ops->restorearrayread) {
4408: (*aij->ops->restorearrayread)(A,array);
4409: } else {
4410: *array = NULL;
4411: }
4412: return 0;
4413: }
4415: /*@C
4416: MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored
4418: Not Collective
4420: Input Parameter:
4421: . mat - a MATSEQAIJ matrix
4423: Output Parameter:
4424: . array - pointer to the data
4426: Level: intermediate
4428: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead()
4429: @*/
4430: PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array)
4431: {
4432: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4434: if (aij->ops->getarraywrite) {
4435: (*aij->ops->getarraywrite)(A,array);
4436: } else {
4437: *array = aij->a;
4438: }
4439: MatSeqAIJInvalidateDiagonal(A);
4440: PetscObjectStateIncrease((PetscObject)A);
4441: return 0;
4442: }
4444: /*@C
4445: MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead
4447: Not Collective
4449: Input Parameter:
4450: . mat - a MATSEQAIJ matrix
4452: Output Parameter:
4453: . array - pointer to the data
4455: Level: intermediate
4457: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4458: @*/
4459: PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array)
4460: {
4461: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4463: if (aij->ops->restorearraywrite) {
4464: (*aij->ops->restorearraywrite)(A,array);
4465: } else {
4466: *array = NULL;
4467: }
4468: return 0;
4469: }
4471: /*@C
4472: MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix
4474: Not Collective
4476: Input Parameter:
4477: . mat - a matrix of type MATSEQAIJ or its subclasses
4479: Output Parameters:
4480: + i - row map array of the matrix
4481: . j - column index array of the matrix
4482: . a - data array of the matrix
4483: - memtype - memory type of the arrays
4485: Notes:
4486: Any of the output parameters can be NULL, in which case the corresponding value is not returned.
4487: If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host.
4489: One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix.
4490: If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix.
4492: Level: Developer
4494: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4495: @*/
4496: PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype)
4497: {
4498: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
4501: if (aij->ops->getcsrandmemtype) {
4502: (*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype);
4503: } else {
4504: if (i) *i = aij->i;
4505: if (j) *j = aij->j;
4506: if (a) *a = aij->a;
4507: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
4508: }
4509: return 0;
4510: }
4512: /*@C
4513: MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row
4515: Not Collective
4517: Input Parameter:
4518: . mat - a MATSEQAIJ matrix
4520: Output Parameter:
4521: . nz - the maximum number of nonzeros in any row
4523: Level: intermediate
4525: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4526: @*/
4527: PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz)
4528: {
4529: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4531: *nz = aij->rmax;
4532: return 0;
4533: }
4535: PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[])
4536: {
4537: MPI_Comm comm;
4538: PetscInt *i,*j;
4539: PetscInt M,N,row;
4540: PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */
4541: PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */
4542: PetscInt *Aj;
4543: PetscScalar *Aa;
4544: Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data);
4545: MatType rtype;
4546: PetscCount *perm,*jmap;
4548: MatResetPreallocationCOO_SeqAIJ(mat);
4549: PetscObjectGetComm((PetscObject)mat,&comm);
4550: MatGetSize(mat,&M,&N);
4551: PetscMalloc2(coo_n,&i,coo_n,&j);
4552: PetscArraycpy(i,coo_i,coo_n); /* Make a copy since we'll modify it */
4553: PetscArraycpy(j,coo_j,coo_n);
4554: PetscMalloc1(coo_n,&perm);
4555: for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */
4556: if (j[k] < 0) i[k] = -1;
4557: perm[k] = k;
4558: }
4560: /* Sort by row */
4561: PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm);
4562: for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */
4563: nneg = k;
4564: PetscMalloc1(coo_n-nneg+1,&jmap); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */
4565: nnz = 0; /* Total number of unique nonzeros to be counted */
4566: jmap++; /* Inc jmap by 1 for convinience */
4568: PetscCalloc1(M+1,&Ai); /* CSR of A */
4569: PetscMalloc1(coo_n-nneg,&Aj); /* We have at most coo_n-nneg unique nonzeros */
4571: /* In each row, sort by column, then unique column indices to get row length */
4572: Ai++; /* Inc by 1 for convinience */
4573: q = 0; /* q-th unique nonzero, with q starting from 0 */
4574: while (k<coo_n) {
4575: row = i[k];
4576: start = k; /* [start,end) indices for this row */
4577: while (k<coo_n && i[k] == row) k++;
4578: end = k;
4579: PetscSortIntWithCountArray(end-start,j+start,perm+start);
4580: /* Find number of unique col entries in this row */
4581: Aj[q] = j[start]; /* Log the first nonzero in this row */
4582: jmap[q] = 1; /* Number of repeats of this nozero entry */
4583: Ai[row] = 1;
4584: nnz++;
4586: for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */
4587: if (j[p] != j[p-1]) { /* Meet a new nonzero */
4588: q++;
4589: jmap[q] = 1;
4590: Aj[q] = j[p];
4591: Ai[row]++;
4592: nnz++;
4593: } else {
4594: jmap[q]++;
4595: }
4596: }
4597: q++; /* Move to next row and thus next unique nonzero */
4598: }
4599: PetscFree2(i,j);
4601: Ai--; /* Back to the beginning of Ai[] */
4602: for (k=0; k<M; k++) Ai[k+1] += Ai[k];
4603: jmap--; /* Back to the beginning of jmap[] */
4604: jmap[0] = 0;
4605: for (k=0; k<nnz; k++) jmap[k+1] += jmap[k];
4606: if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */
4607: PetscCount *jmap_new;
4608: PetscInt *Aj_new;
4610: PetscMalloc1(nnz+1,&jmap_new);
4611: PetscArraycpy(jmap_new,jmap,nnz+1);
4612: PetscFree(jmap);
4613: jmap = jmap_new;
4615: PetscMalloc1(nnz,&Aj_new);
4616: PetscArraycpy(Aj_new,Aj,nnz);
4617: PetscFree(Aj);
4618: Aj = Aj_new;
4619: }
4621: if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */
4622: PetscCount *perm_new;
4624: PetscMalloc1(coo_n-nneg,&perm_new);
4625: PetscArraycpy(perm_new,perm+nneg,coo_n-nneg);
4626: PetscFree(perm);
4627: perm = perm_new;
4628: }
4630: MatGetRootType_Private(mat,&rtype);
4631: PetscCalloc1(nnz,&Aa); /* Zero the matrix */
4632: MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat);
4634: seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */
4635: seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */
4636: /* Record COO fields */
4637: seqaij->coo_n = coo_n;
4638: seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */
4639: seqaij->jmap = jmap; /* of length nnz+1 */
4640: seqaij->perm = perm;
4641: return 0;
4642: }
4644: static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode)
4645: {
4646: Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data;
4647: PetscCount i,j,Annz = aseq->nz;
4648: PetscCount *perm = aseq->perm,*jmap = aseq->jmap;
4649: PetscScalar *Aa;
4651: MatSeqAIJGetArray(A,&Aa);
4652: for (i=0; i<Annz; i++) {
4653: PetscScalar sum = 0.0;
4654: for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]];
4655: Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum;
4656: }
4657: MatSeqAIJRestoreArray(A,&Aa);
4658: return 0;
4659: }
4661: #if defined(PETSC_HAVE_CUDA)
4662: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*);
4663: #endif
4664: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4665: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*);
4666: #endif
4668: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B)
4669: {
4670: Mat_SeqAIJ *b;
4671: PetscMPIInt size;
4673: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
4676: PetscNewLog(B,&b);
4678: B->data = (void*)b;
4680: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
4681: if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull;
4683: b->row = NULL;
4684: b->col = NULL;
4685: b->icol = NULL;
4686: b->reallocs = 0;
4687: b->ignorezeroentries = PETSC_FALSE;
4688: b->roworiented = PETSC_TRUE;
4689: b->nonew = 0;
4690: b->diag = NULL;
4691: b->solve_work = NULL;
4692: B->spptr = NULL;
4693: b->saved_values = NULL;
4694: b->idiag = NULL;
4695: b->mdiag = NULL;
4696: b->ssor_work = NULL;
4697: b->omega = 1.0;
4698: b->fshift = 0.0;
4699: b->idiagvalid = PETSC_FALSE;
4700: b->ibdiagvalid = PETSC_FALSE;
4701: b->keepnonzeropattern = PETSC_FALSE;
4703: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4705: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4706: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);
4707: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);
4708: #endif
4710: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);
4711: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);
4712: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);
4713: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);
4714: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);
4715: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);
4716: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);
4717: #if defined(PETSC_HAVE_MKL_SPARSE)
4718: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);
4719: #endif
4720: #if defined(PETSC_HAVE_CUDA)
4721: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);
4722: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4723: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ);
4724: #endif
4725: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4726: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos);
4727: #endif
4728: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);
4729: #if defined(PETSC_HAVE_ELEMENTAL)
4730: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);
4731: #endif
4732: #if defined(PETSC_HAVE_SCALAPACK)
4733: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK);
4734: #endif
4735: #if defined(PETSC_HAVE_HYPRE)
4736: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);
4737: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ);
4738: #endif
4739: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);
4740: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);
4741: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);
4742: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);
4743: PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);
4744: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);
4745: PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);
4746: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);
4747: PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);
4748: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ);
4749: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ);
4750: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4751: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ);
4752: PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ);
4753: PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ);
4754: MatCreate_SeqAIJ_Inode(B);
4755: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4756: MatSeqAIJSetTypeFromOptions(B); /* this allows changing the matrix subtype to say MATSEQAIJPERM */
4757: return 0;
4758: }
4760: /*
4761: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
4762: */
4763: PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
4764: {
4765: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data;
4766: PetscInt m = A->rmap->n,i;
4770: C->factortype = A->factortype;
4771: c->row = NULL;
4772: c->col = NULL;
4773: c->icol = NULL;
4774: c->reallocs = 0;
4776: C->assembled = A->assembled;
4777: C->preallocated = A->preallocated;
4779: if (A->preallocated) {
4780: PetscLayoutReference(A->rmap,&C->rmap);
4781: PetscLayoutReference(A->cmap,&C->cmap);
4783: PetscMalloc1(m,&c->imax);
4784: PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt));
4785: PetscMalloc1(m,&c->ilen);
4786: PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt));
4787: PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));
4789: /* allocate the matrix space */
4790: if (mallocmatspace) {
4791: PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);
4792: PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));
4794: c->singlemalloc = PETSC_TRUE;
4796: PetscArraycpy(c->i,a->i,m+1);
4797: if (m > 0) {
4798: PetscArraycpy(c->j,a->j,a->i[m]);
4799: if (cpvalues == MAT_COPY_VALUES) {
4800: const PetscScalar *aa;
4802: MatSeqAIJGetArrayRead(A,&aa);
4803: PetscArraycpy(c->a,aa,a->i[m]);
4804: MatSeqAIJGetArrayRead(A,&aa);
4805: } else {
4806: PetscArrayzero(c->a,a->i[m]);
4807: }
4808: }
4809: }
4811: c->ignorezeroentries = a->ignorezeroentries;
4812: c->roworiented = a->roworiented;
4813: c->nonew = a->nonew;
4814: if (a->diag) {
4815: PetscMalloc1(m+1,&c->diag);
4816: PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt));
4817: PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));
4818: } else c->diag = NULL;
4820: c->solve_work = NULL;
4821: c->saved_values = NULL;
4822: c->idiag = NULL;
4823: c->ssor_work = NULL;
4824: c->keepnonzeropattern = a->keepnonzeropattern;
4825: c->free_a = PETSC_TRUE;
4826: c->free_ij = PETSC_TRUE;
4828: c->rmax = a->rmax;
4829: c->nz = a->nz;
4830: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
4832: c->compressedrow.use = a->compressedrow.use;
4833: c->compressedrow.nrows = a->compressedrow.nrows;
4834: if (a->compressedrow.use) {
4835: i = a->compressedrow.nrows;
4836: PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);
4837: PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);
4838: PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);
4839: } else {
4840: c->compressedrow.use = PETSC_FALSE;
4841: c->compressedrow.i = NULL;
4842: c->compressedrow.rindex = NULL;
4843: }
4844: c->nonzerorowcnt = a->nonzerorowcnt;
4845: C->nonzerostate = A->nonzerostate;
4847: MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);
4848: }
4849: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
4850: return 0;
4851: }
4853: PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
4854: {
4855: MatCreate(PetscObjectComm((PetscObject)A),B);
4856: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
4857: if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) {
4858: MatSetBlockSizesFromMats(*B,A,A);
4859: }
4860: MatSetType(*B,((PetscObject)A)->type_name);
4861: MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);
4862: return 0;
4863: }
4865: PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer)
4866: {
4867: PetscBool isbinary, ishdf5;
4871: /* force binary viewer to load .info file if it has not yet done so */
4872: PetscViewerSetUp(viewer);
4873: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
4874: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
4875: if (isbinary) {
4876: MatLoad_SeqAIJ_Binary(newMat,viewer);
4877: } else if (ishdf5) {
4878: #if defined(PETSC_HAVE_HDF5)
4879: MatLoad_AIJ_HDF5(newMat,viewer);
4880: #else
4881: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
4882: #endif
4883: } else {
4884: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name);
4885: }
4886: return 0;
4887: }
4889: PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer)
4890: {
4891: Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data;
4892: PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i;
4894: PetscViewerSetUp(viewer);
4896: /* read in matrix header */
4897: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
4899: M = header[1]; N = header[2]; nz = header[3];
4904: /* set block sizes from the viewer's .info file */
4905: MatLoad_Binary_BlockSizes(mat,viewer);
4906: /* set local and global sizes if not set already */
4907: if (mat->rmap->n < 0) mat->rmap->n = M;
4908: if (mat->cmap->n < 0) mat->cmap->n = N;
4909: if (mat->rmap->N < 0) mat->rmap->N = M;
4910: if (mat->cmap->N < 0) mat->cmap->N = N;
4911: PetscLayoutSetUp(mat->rmap);
4912: PetscLayoutSetUp(mat->cmap);
4914: /* check if the matrix sizes are correct */
4915: MatGetSize(mat,&rows,&cols);
4918: /* read in row lengths */
4919: PetscMalloc1(M,&rowlens);
4920: PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT);
4921: /* check if sum(rowlens) is same as nz */
4922: sum = 0; for (i=0; i<M; i++) sum += rowlens[i];
4924: /* preallocate and check sizes */
4925: MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens);
4926: MatGetSize(mat,&rows,&cols);
4928: /* store row lengths */
4929: PetscArraycpy(a->ilen,rowlens,M);
4930: PetscFree(rowlens);
4932: /* fill in "i" row pointers */
4933: a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i];
4934: /* read in "j" column indices */
4935: PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT);
4936: /* read in "a" nonzero values */
4937: PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR);
4939: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
4940: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
4941: return 0;
4942: }
4944: PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg)
4945: {
4946: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data;
4947: const PetscScalar *aa,*ba;
4948: #if defined(PETSC_USE_COMPLEX)
4949: PetscInt k;
4950: #endif
4952: /* If the matrix dimensions are not equal,or no of nonzeros */
4953: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) {
4954: *flg = PETSC_FALSE;
4955: return 0;
4956: }
4958: /* if the a->i are the same */
4959: PetscArraycmp(a->i,b->i,A->rmap->n+1,flg);
4960: if (!*flg) return 0;
4962: /* if a->j are the same */
4963: PetscArraycmp(a->j,b->j,a->nz,flg);
4964: if (!*flg) return 0;
4966: MatSeqAIJGetArrayRead(A,&aa);
4967: MatSeqAIJGetArrayRead(B,&ba);
4968: /* if a->a are the same */
4969: #if defined(PETSC_USE_COMPLEX)
4970: for (k=0; k<a->nz; k++) {
4971: if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) {
4972: *flg = PETSC_FALSE;
4973: return 0;
4974: }
4975: }
4976: #else
4977: PetscArraycmp(aa,ba,a->nz,flg);
4978: #endif
4979: MatSeqAIJRestoreArrayRead(A,&aa);
4980: MatSeqAIJRestoreArrayRead(B,&ba);
4981: return 0;
4982: }
4984: /*@
4985: MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format)
4986: provided by the user.
4988: Collective
4990: Input Parameters:
4991: + comm - must be an MPI communicator of size 1
4992: . m - number of rows
4993: . n - number of columns
4994: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix
4995: . j - column indices
4996: - a - matrix values
4998: Output Parameter:
4999: . mat - the matrix
5001: Level: intermediate
5003: Notes:
5004: The i, j, and a arrays are not copied by this routine, the user must free these arrays
5005: once the matrix is destroyed and not before
5007: You cannot set new nonzero locations into this matrix, that will generate an error.
5009: The i and j indices are 0 based
5011: The format which is used for the sparse matrix input, is equivalent to a
5012: row-major ordering.. i.e for the following matrix, the input data expected is
5013: as shown
5015: $ 1 0 0
5016: $ 2 0 3
5017: $ 4 5 6
5018: $
5019: $ i = {0,1,3,6} [size = nrow+1 = 3+1]
5020: $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row
5021: $ v = {1,2,3,4,5,6} [size = 6]
5023: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5025: @*/
5026: PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
5027: {
5028: PetscInt ii;
5029: Mat_SeqAIJ *aij;
5030: PetscInt jj;
5033: MatCreate(comm,mat);
5034: MatSetSizes(*mat,m,n,m,n);
5035: /* MatSetBlockSizes(*mat,,); */
5036: MatSetType(*mat,MATSEQAIJ);
5037: MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL);
5038: aij = (Mat_SeqAIJ*)(*mat)->data;
5039: PetscMalloc1(m,&aij->imax);
5040: PetscMalloc1(m,&aij->ilen);
5042: aij->i = i;
5043: aij->j = j;
5044: aij->a = a;
5045: aij->singlemalloc = PETSC_FALSE;
5046: aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
5047: aij->free_a = PETSC_FALSE;
5048: aij->free_ij = PETSC_FALSE;
5050: for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) {
5051: aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
5052: if (PetscDefined(USE_DEBUG)) {
5054: for (jj=i[ii]+1; jj<i[ii+1]; jj++) {
5057: }
5058: }
5059: }
5060: if (PetscDefined(USE_DEBUG)) {
5061: for (ii=0; ii<aij->i[m]; ii++) {
5064: }
5065: }
5067: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5068: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5069: return 0;
5070: }
5072: /*@C
5073: MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format)
5074: provided by the user.
5076: Collective
5078: Input Parameters:
5079: + comm - must be an MPI communicator of size 1
5080: . m - number of rows
5081: . n - number of columns
5082: . i - row indices
5083: . j - column indices
5084: . a - matrix values
5085: . nz - number of nonzeros
5086: - idx - 0 or 1 based
5088: Output Parameter:
5089: . mat - the matrix
5091: Level: intermediate
5093: Notes:
5094: The i and j indices are 0 based. The format which is used for the sparse matrix input, is equivalent to a row-major ordering. i.e for the following matrix,
5095: the input data expected is as shown
5096: .vb
5097: 1 0 0
5098: 2 0 3
5099: 4 5 6
5101: i = {0,1,1,2,2,2}
5102: j = {0,0,2,0,1,2}
5103: v = {1,2,3,4,5,6}
5104: .ve
5106: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5108: @*/
5109: PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx)
5110: {
5111: PetscInt ii, *nnz, one = 1,row,col;
5113: PetscCalloc1(m,&nnz);
5114: for (ii = 0; ii < nz; ii++) {
5115: nnz[i[ii] - !!idx] += 1;
5116: }
5117: MatCreate(comm,mat);
5118: MatSetSizes(*mat,m,n,m,n);
5119: MatSetType(*mat,MATSEQAIJ);
5120: MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);
5121: for (ii = 0; ii < nz; ii++) {
5122: if (idx) {
5123: row = i[ii] - 1;
5124: col = j[ii] - 1;
5125: } else {
5126: row = i[ii];
5127: col = j[ii];
5128: }
5129: MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);
5130: }
5131: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5132: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5133: PetscFree(nnz);
5134: return 0;
5135: }
5137: PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A)
5138: {
5139: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
5141: a->idiagvalid = PETSC_FALSE;
5142: a->ibdiagvalid = PETSC_FALSE;
5144: MatSeqAIJInvalidateDiagonal_Inode(A);
5145: return 0;
5146: }
5148: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
5149: {
5150: PetscMPIInt size;
5152: MPI_Comm_size(comm,&size);
5153: if (size == 1) {
5154: if (scall == MAT_INITIAL_MATRIX) {
5155: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
5156: } else {
5157: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
5158: }
5159: } else {
5160: MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);
5161: }
5162: return 0;
5163: }
5165: /*
5166: Permute A into C's *local* index space using rowemb,colemb.
5167: The embedding are supposed to be injections and the above implies that the range of rowemb is a subset
5168: of [0,m), colemb is in [0,n).
5169: If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A.
5170: */
5171: PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B)
5172: {
5173: /* If making this function public, change the error returned in this function away from _PLIB. */
5174: Mat_SeqAIJ *Baij;
5175: PetscBool seqaij;
5176: PetscInt m,n,*nz,i,j,count;
5177: PetscScalar v;
5178: const PetscInt *rowindices,*colindices;
5180: if (!B) return 0;
5181: /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */
5182: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
5184: if (rowemb) {
5185: ISGetLocalSize(rowemb,&m);
5187: } else {
5189: }
5190: if (colemb) {
5191: ISGetLocalSize(colemb,&n);
5193: } else {
5195: }
5197: Baij = (Mat_SeqAIJ*)(B->data);
5198: if (pattern == DIFFERENT_NONZERO_PATTERN) {
5199: PetscMalloc1(B->rmap->n,&nz);
5200: for (i=0; i<B->rmap->n; i++) {
5201: nz[i] = Baij->i[i+1] - Baij->i[i];
5202: }
5203: MatSeqAIJSetPreallocation(C,0,nz);
5204: PetscFree(nz);
5205: }
5206: if (pattern == SUBSET_NONZERO_PATTERN) {
5207: MatZeroEntries(C);
5208: }
5209: count = 0;
5210: rowindices = NULL;
5211: colindices = NULL;
5212: if (rowemb) {
5213: ISGetIndices(rowemb,&rowindices);
5214: }
5215: if (colemb) {
5216: ISGetIndices(colemb,&colindices);
5217: }
5218: for (i=0; i<B->rmap->n; i++) {
5219: PetscInt row;
5220: row = i;
5221: if (rowindices) row = rowindices[i];
5222: for (j=Baij->i[i]; j<Baij->i[i+1]; j++) {
5223: PetscInt col;
5224: col = Baij->j[count];
5225: if (colindices) col = colindices[col];
5226: v = Baij->a[count];
5227: MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);
5228: ++count;
5229: }
5230: }
5231: /* FIXME: set C's nonzerostate correctly. */
5232: /* Assembly for C is necessary. */
5233: C->preallocated = PETSC_TRUE;
5234: C->assembled = PETSC_TRUE;
5235: C->was_assembled = PETSC_FALSE;
5236: return 0;
5237: }
5239: PetscFunctionList MatSeqAIJList = NULL;
5241: /*@C
5242: MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype
5244: Collective on Mat
5246: Input Parameters:
5247: + mat - the matrix object
5248: - matype - matrix type
5250: Options Database Key:
5251: . -mat_seqai_type <method> - for example seqaijcrl
5253: Level: intermediate
5255: .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat
5256: @*/
5257: PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype)
5258: {
5259: PetscBool sametype;
5260: PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*);
5263: PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);
5264: if (sametype) return 0;
5266: PetscFunctionListFind(MatSeqAIJList,matype,&r);
5268: (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);
5269: return 0;
5270: }
5272: /*@C
5273: MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices
5275: Not Collective
5277: Input Parameters:
5278: + name - name of a new user-defined matrix type, for example MATSEQAIJCRL
5279: - function - routine to convert to subtype
5281: Notes:
5282: MatSeqAIJRegister() may be called multiple times to add several user-defined solvers.
5284: Then, your matrix can be chosen with the procedural interface at runtime via the option
5285: $ -mat_seqaij_type my_mat
5287: Level: advanced
5289: .seealso: MatSeqAIJRegisterAll()
5291: Level: advanced
5292: @*/
5293: PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *))
5294: {
5295: MatInitializePackage();
5296: PetscFunctionListAdd(&MatSeqAIJList,sname,function);
5297: return 0;
5298: }
5300: PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE;
5302: /*@C
5303: MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ
5305: Not Collective
5307: Level: advanced
5309: .seealso: MatRegisterAll(), MatSeqAIJRegister()
5310: @*/
5311: PetscErrorCode MatSeqAIJRegisterAll(void)
5312: {
5313: if (MatSeqAIJRegisterAllCalled) return 0;
5314: MatSeqAIJRegisterAllCalled = PETSC_TRUE;
5316: MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);
5317: MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);
5318: MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);
5319: #if defined(PETSC_HAVE_MKL_SPARSE)
5320: MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);
5321: #endif
5322: #if defined(PETSC_HAVE_CUDA)
5323: MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE);
5324: #endif
5325: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
5326: MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos);
5327: #endif
5328: #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA)
5329: MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);
5330: #endif
5331: return 0;
5332: }
5334: /*
5335: Special version for direct calls from Fortran
5336: */
5337: #include <petsc/private/fortranimpl.h>
5338: #if defined(PETSC_HAVE_FORTRAN_CAPS)
5339: #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ
5340: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
5341: #define matsetvaluesseqaij_ matsetvaluesseqaij
5342: #endif
5344: /* Change these macros so can be used in void function */
5346: /* Change these macros so can be used in void function */
5347: /* Identical to PetscCallVoid, except it assigns to *_ierr */
5348: #undef PetscCall
5349: #define PetscCall(...) do { \
5350: PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \
5351: if (PetscUnlikely(ierr_msv_mpiaij)) { \
5352: *_PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \
5353: return; \
5354: } \
5355: } while (0)
5357: #undef SETERRQ
5358: #define SETERRQ(comm,ierr,...) do { \
5359: *_PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \
5360: return; \
5361: } while (0)
5363: PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr)
5364: {
5365: Mat A = *AA;
5366: PetscInt m = *mm, n = *nn;
5367: InsertMode is = *isis;
5368: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
5369: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
5370: PetscInt *imax,*ai,*ailen;
5371: PetscInt *aj,nonew = a->nonew,lastcol = -1;
5372: MatScalar *ap,value,*aa;
5373: PetscBool ignorezeroentries = a->ignorezeroentries;
5374: PetscBool roworiented = a->roworiented;
5376: MatCheckPreallocated(A,1);
5377: imax = a->imax;
5378: ai = a->i;
5379: ailen = a->ilen;
5380: aj = a->j;
5381: aa = a->a;
5383: for (k=0; k<m; k++) { /* loop over added rows */
5384: row = im[k];
5385: if (row < 0) continue;
5387: rp = aj + ai[row]; ap = aa + ai[row];
5388: rmax = imax[row]; nrow = ailen[row];
5389: low = 0;
5390: high = nrow;
5391: for (l=0; l<n; l++) { /* loop over added columns */
5392: if (in[l] < 0) continue;
5394: col = in[l];
5395: if (roworiented) value = v[l + k*n];
5396: else value = v[k + l*m];
5398: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
5400: if (col <= lastcol) low = 0;
5401: else high = nrow;
5402: lastcol = col;
5403: while (high-low > 5) {
5404: t = (low+high)/2;
5405: if (rp[t] > col) high = t;
5406: else low = t;
5407: }
5408: for (i=low; i<high; i++) {
5409: if (rp[i] > col) break;
5410: if (rp[i] == col) {
5411: if (is == ADD_VALUES) ap[i] += value;
5412: else ap[i] = value;
5413: goto noinsert;
5414: }
5415: }
5416: if (value == 0.0 && ignorezeroentries) goto noinsert;
5417: if (nonew == 1) goto noinsert;
5419: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
5420: N = nrow++ - 1; a->nz++; high++;
5421: /* shift up all the later entries in this row */
5422: for (ii=N; ii>=i; ii--) {
5423: rp[ii+1] = rp[ii];
5424: ap[ii+1] = ap[ii];
5425: }
5426: rp[i] = col;
5427: ap[i] = value;
5428: A->nonzerostate++;
5429: noinsert:;
5430: low = i + 1;
5431: }
5432: ailen[row] = nrow;
5433: }
5434: return;
5435: }
5436: /* Undefining these here since they were redefined from their original definition above! No
5437: * other PETSc functions should be defined past this point, as it is impossible to recover the
5438: * original definitions */
5439: #undef PetscCall
5440: #undef SETERRQ