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 (i < A->cmap->n && diag[i] >= ii[i+1]) { /* 'out of range' rows never have diagonals */
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<PetscMin(A->rmap->n,A->cmap->n); i++) {
1712: a->imax[i] += mdiag[i];
1713: }
1714: MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax);
1716: /* copy old values into new matrix data structure */
1717: for (i=0; i<A->rmap->n; i++) {
1718: MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES);
1719: if (i < A->cmap->n) {
1720: MatSetValue(A,i,i,v,ADD_VALUES);
1721: }
1722: }
1723: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1724: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1725: if (singlemalloc) {
1726: PetscFree3(olda,oldj,oldi);
1727: } else {
1728: if (free_a) PetscFree(olda);
1729: if (free_ij) PetscFree(oldj);
1730: if (free_ij) PetscFree(oldi);
1731: }
1732: }
1733: PetscFree(mdiag);
1734: a->diagonaldense = PETSC_TRUE;
1735: return 0;
1736: }
1738: /*
1739: Checks for missing diagonals
1740: */
1741: PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d)
1742: {
1743: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1744: PetscInt *diag,*ii = a->i,i;
1746: *missing = PETSC_FALSE;
1747: if (A->rmap->n > 0 && !ii) {
1748: *missing = PETSC_TRUE;
1749: if (d) *d = 0;
1750: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
1751: } else {
1752: PetscInt n;
1753: n = PetscMin(A->rmap->n, A->cmap->n);
1754: diag = a->diag;
1755: for (i=0; i<n; i++) {
1756: if (diag[i] >= ii[i+1]) {
1757: *missing = PETSC_TRUE;
1758: if (d) *d = i;
1759: PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i);
1760: break;
1761: }
1762: }
1763: }
1764: return 0;
1765: }
1767: #include <petscblaslapack.h>
1768: #include <petsc/private/kernels/blockinvert.h>
1770: /*
1771: Note that values is allocated externally by the PC and then passed into this routine
1772: */
1773: PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag)
1774: {
1775: PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots;
1776: PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
1777: const PetscReal shift = 0.0;
1778: PetscInt ipvt[5];
1779: PetscScalar work[25],*v_work;
1781: allowzeropivot = PetscNot(A->erroriffailure);
1782: for (i=0; i<nblocks; i++) ncnt += bsizes[i];
1784: for (i=0; i<nblocks; i++) {
1785: bsizemax = PetscMax(bsizemax,bsizes[i]);
1786: }
1787: PetscMalloc1(bsizemax,&indx);
1788: if (bsizemax > 7) {
1789: PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots);
1790: }
1791: ncnt = 0;
1792: for (i=0; i<nblocks; i++) {
1793: for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j;
1794: MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag);
1795: switch (bsizes[i]) {
1796: case 1:
1797: *diag = 1.0/(*diag);
1798: break;
1799: case 2:
1800: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
1801: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1802: PetscKernel_A_gets_transpose_A_2(diag);
1803: break;
1804: case 3:
1805: PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
1806: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1807: PetscKernel_A_gets_transpose_A_3(diag);
1808: break;
1809: case 4:
1810: PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
1811: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1812: PetscKernel_A_gets_transpose_A_4(diag);
1813: break;
1814: case 5:
1815: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
1816: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1817: PetscKernel_A_gets_transpose_A_5(diag);
1818: break;
1819: case 6:
1820: PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
1821: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1822: PetscKernel_A_gets_transpose_A_6(diag);
1823: break;
1824: case 7:
1825: PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
1826: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1827: PetscKernel_A_gets_transpose_A_7(diag);
1828: break;
1829: default:
1830: PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
1831: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1832: PetscKernel_A_gets_transpose_A_N(diag,bsizes[i]);
1833: }
1834: ncnt += bsizes[i];
1835: diag += bsizes[i]*bsizes[i];
1836: }
1837: if (bsizemax > 7) {
1838: PetscFree2(v_work,v_pivots);
1839: }
1840: PetscFree(indx);
1841: return 0;
1842: }
1844: /*
1845: Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways
1846: */
1847: PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift)
1848: {
1849: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
1850: PetscInt i,*diag,m = A->rmap->n;
1851: const MatScalar *v;
1852: PetscScalar *idiag,*mdiag;
1854: if (a->idiagvalid) return 0;
1855: MatMarkDiagonal_SeqAIJ(A);
1856: diag = a->diag;
1857: if (!a->idiag) {
1858: PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);
1859: PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar));
1860: }
1862: mdiag = a->mdiag;
1863: idiag = a->idiag;
1864: MatSeqAIJGetArrayRead(A,&v);
1865: if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) {
1866: for (i=0; i<m; i++) {
1867: mdiag[i] = v[diag[i]];
1868: if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */
1869: if (PetscRealPart(fshift)) {
1870: PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i);
1871: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1872: A->factorerror_zeropivot_value = 0.0;
1873: A->factorerror_zeropivot_row = i;
1874: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i);
1875: }
1876: idiag[i] = 1.0/v[diag[i]];
1877: }
1878: PetscLogFlops(m);
1879: } else {
1880: for (i=0; i<m; i++) {
1881: mdiag[i] = v[diag[i]];
1882: idiag[i] = omega/(fshift + v[diag[i]]);
1883: }
1884: PetscLogFlops(2.0*m);
1885: }
1886: a->idiagvalid = PETSC_TRUE;
1887: MatSeqAIJRestoreArrayRead(A,&v);
1888: return 0;
1889: }
1891: #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h>
1892: PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1893: {
1894: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1895: PetscScalar *x,d,sum,*t,scale;
1896: const MatScalar *v,*idiag=NULL,*mdiag,*aa;
1897: const PetscScalar *b, *bs,*xb, *ts;
1898: PetscInt n,m = A->rmap->n,i;
1899: const PetscInt *idx,*diag;
1901: if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) {
1902: MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx);
1903: return 0;
1904: }
1905: its = its*lits;
1907: if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */
1908: if (!a->idiagvalid) MatInvertDiagonal_SeqAIJ(A,omega,fshift);
1909: a->fshift = fshift;
1910: a->omega = omega;
1912: diag = a->diag;
1913: t = a->ssor_work;
1914: idiag = a->idiag;
1915: mdiag = a->mdiag;
1917: MatSeqAIJGetArrayRead(A,&aa);
1918: VecGetArray(xx,&x);
1919: VecGetArrayRead(bb,&b);
1920: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
1921: if (flag == SOR_APPLY_UPPER) {
1922: /* apply (U + D/omega) to the vector */
1923: bs = b;
1924: for (i=0; i<m; i++) {
1925: d = fshift + mdiag[i];
1926: n = a->i[i+1] - diag[i] - 1;
1927: idx = a->j + diag[i] + 1;
1928: v = aa + diag[i] + 1;
1929: sum = b[i]*d/omega;
1930: PetscSparseDensePlusDot(sum,bs,v,idx,n);
1931: x[i] = sum;
1932: }
1933: VecRestoreArray(xx,&x);
1934: VecRestoreArrayRead(bb,&b);
1935: MatSeqAIJRestoreArrayRead(A,&aa);
1936: PetscLogFlops(a->nz);
1937: return 0;
1938: }
1941: else if (flag & SOR_EISENSTAT) {
1942: /* Let A = L + U + D; where L is lower triangular,
1943: U is upper triangular, E = D/omega; This routine applies
1945: (L + E)^{-1} A (U + E)^{-1}
1947: to a vector efficiently using Eisenstat's trick.
1948: */
1949: scale = (2.0/omega) - 1.0;
1951: /* x = (E + U)^{-1} b */
1952: for (i=m-1; i>=0; i--) {
1953: n = a->i[i+1] - diag[i] - 1;
1954: idx = a->j + diag[i] + 1;
1955: v = aa + diag[i] + 1;
1956: sum = b[i];
1957: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1958: x[i] = sum*idiag[i];
1959: }
1961: /* t = b - (2*E - D)x */
1962: v = aa;
1963: for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i];
1965: /* t = (E + L)^{-1}t */
1966: ts = t;
1967: diag = a->diag;
1968: for (i=0; i<m; i++) {
1969: n = diag[i] - a->i[i];
1970: idx = a->j + a->i[i];
1971: v = aa + a->i[i];
1972: sum = t[i];
1973: PetscSparseDenseMinusDot(sum,ts,v,idx,n);
1974: t[i] = sum*idiag[i];
1975: /* x = x + t */
1976: x[i] += t[i];
1977: }
1979: PetscLogFlops(6.0*m-1 + 2.0*a->nz);
1980: VecRestoreArray(xx,&x);
1981: VecRestoreArrayRead(bb,&b);
1982: return 0;
1983: }
1984: if (flag & SOR_ZERO_INITIAL_GUESS) {
1985: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1986: for (i=0; i<m; i++) {
1987: n = diag[i] - a->i[i];
1988: idx = a->j + a->i[i];
1989: v = aa + a->i[i];
1990: sum = b[i];
1991: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1992: t[i] = sum;
1993: x[i] = sum*idiag[i];
1994: }
1995: xb = t;
1996: PetscLogFlops(a->nz);
1997: } else xb = b;
1998: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1999: for (i=m-1; i>=0; i--) {
2000: n = a->i[i+1] - diag[i] - 1;
2001: idx = a->j + diag[i] + 1;
2002: v = aa + diag[i] + 1;
2003: sum = xb[i];
2004: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2005: if (xb == b) {
2006: x[i] = sum*idiag[i];
2007: } else {
2008: x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2009: }
2010: }
2011: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
2012: }
2013: its--;
2014: }
2015: while (its--) {
2016: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
2017: for (i=0; i<m; i++) {
2018: /* lower */
2019: n = diag[i] - a->i[i];
2020: idx = a->j + a->i[i];
2021: v = aa + a->i[i];
2022: sum = b[i];
2023: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2024: t[i] = sum; /* save application of the lower-triangular part */
2025: /* upper */
2026: n = a->i[i+1] - diag[i] - 1;
2027: idx = a->j + diag[i] + 1;
2028: v = aa + diag[i] + 1;
2029: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2030: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2031: }
2032: xb = t;
2033: PetscLogFlops(2.0*a->nz);
2034: } else xb = b;
2035: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
2036: for (i=m-1; i>=0; i--) {
2037: sum = xb[i];
2038: if (xb == b) {
2039: /* whole matrix (no checkpointing available) */
2040: n = a->i[i+1] - a->i[i];
2041: idx = a->j + a->i[i];
2042: v = aa + a->i[i];
2043: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2044: x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i];
2045: } else { /* lower-triangular part has been saved, so only apply upper-triangular */
2046: n = a->i[i+1] - diag[i] - 1;
2047: idx = a->j + diag[i] + 1;
2048: v = aa + diag[i] + 1;
2049: PetscSparseDenseMinusDot(sum,x,v,idx,n);
2050: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
2051: }
2052: }
2053: if (xb == b) {
2054: PetscLogFlops(2.0*a->nz);
2055: } else {
2056: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
2057: }
2058: }
2059: }
2060: MatSeqAIJRestoreArrayRead(A,&aa);
2061: VecRestoreArray(xx,&x);
2062: VecRestoreArrayRead(bb,&b);
2063: return 0;
2064: }
2066: PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info)
2067: {
2068: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2070: info->block_size = 1.0;
2071: info->nz_allocated = a->maxnz;
2072: info->nz_used = a->nz;
2073: info->nz_unneeded = (a->maxnz - a->nz);
2074: info->assemblies = A->num_ass;
2075: info->mallocs = A->info.mallocs;
2076: info->memory = ((PetscObject)A)->mem;
2077: if (A->factortype) {
2078: info->fill_ratio_given = A->info.fill_ratio_given;
2079: info->fill_ratio_needed = A->info.fill_ratio_needed;
2080: info->factor_mallocs = A->info.factor_mallocs;
2081: } else {
2082: info->fill_ratio_given = 0;
2083: info->fill_ratio_needed = 0;
2084: info->factor_mallocs = 0;
2085: }
2086: return 0;
2087: }
2089: PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2090: {
2091: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2092: PetscInt i,m = A->rmap->n - 1;
2093: const PetscScalar *xx;
2094: PetscScalar *bb,*aa;
2095: PetscInt d = 0;
2097: if (x && b) {
2098: VecGetArrayRead(x,&xx);
2099: VecGetArray(b,&bb);
2100: for (i=0; i<N; i++) {
2102: if (rows[i] >= A->cmap->n) continue;
2103: bb[rows[i]] = diag*xx[rows[i]];
2104: }
2105: VecRestoreArrayRead(x,&xx);
2106: VecRestoreArray(b,&bb);
2107: }
2109: MatSeqAIJGetArray(A,&aa);
2110: if (a->keepnonzeropattern) {
2111: for (i=0; i<N; i++) {
2113: PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]]);
2114: }
2115: if (diag != 0.0) {
2116: for (i=0; i<N; i++) {
2117: d = rows[i];
2118: if (rows[i] >= A->cmap->n) continue;
2120: }
2121: for (i=0; i<N; i++) {
2122: if (rows[i] >= A->cmap->n) continue;
2123: aa[a->diag[rows[i]]] = diag;
2124: }
2125: }
2126: } else {
2127: if (diag != 0.0) {
2128: for (i=0; i<N; i++) {
2130: if (a->ilen[rows[i]] > 0) {
2131: if (rows[i] >= A->cmap->n) {
2132: a->ilen[rows[i]] = 0;
2133: } else {
2134: a->ilen[rows[i]] = 1;
2135: aa[a->i[rows[i]]] = diag;
2136: a->j[a->i[rows[i]]] = rows[i];
2137: }
2138: } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */
2139: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
2140: }
2141: }
2142: } else {
2143: for (i=0; i<N; i++) {
2145: a->ilen[rows[i]] = 0;
2146: }
2147: }
2148: A->nonzerostate++;
2149: }
2150: MatSeqAIJRestoreArray(A,&aa);
2151: (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);
2152: return 0;
2153: }
2155: PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2156: {
2157: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2158: PetscInt i,j,m = A->rmap->n - 1,d = 0;
2159: PetscBool missing,*zeroed,vecs = PETSC_FALSE;
2160: const PetscScalar *xx;
2161: PetscScalar *bb,*aa;
2163: if (!N) return 0;
2164: MatSeqAIJGetArray(A,&aa);
2165: if (x && b) {
2166: VecGetArrayRead(x,&xx);
2167: VecGetArray(b,&bb);
2168: vecs = PETSC_TRUE;
2169: }
2170: PetscCalloc1(A->rmap->n,&zeroed);
2171: for (i=0; i<N; i++) {
2173: PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]]);
2175: zeroed[rows[i]] = PETSC_TRUE;
2176: }
2177: for (i=0; i<A->rmap->n; i++) {
2178: if (!zeroed[i]) {
2179: for (j=a->i[i]; j<a->i[i+1]; j++) {
2180: if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) {
2181: if (vecs) bb[i] -= aa[j]*xx[a->j[j]];
2182: aa[j] = 0.0;
2183: }
2184: }
2185: } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i];
2186: }
2187: if (x && b) {
2188: VecRestoreArrayRead(x,&xx);
2189: VecRestoreArray(b,&bb);
2190: }
2191: PetscFree(zeroed);
2192: if (diag != 0.0) {
2193: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
2194: if (missing) {
2195: for (i=0; i<N; i++) {
2196: if (rows[i] >= A->cmap->N) continue;
2198: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
2199: }
2200: } else {
2201: for (i=0; i<N; i++) {
2202: aa[a->diag[rows[i]]] = diag;
2203: }
2204: }
2205: }
2206: MatSeqAIJRestoreArray(A,&aa);
2207: (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);
2208: return 0;
2209: }
2211: PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
2212: {
2213: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2214: const PetscScalar *aa;
2215: PetscInt *itmp;
2217: MatSeqAIJGetArrayRead(A,&aa);
2218: *nz = a->i[row+1] - a->i[row];
2219: if (v) *v = (PetscScalar*)(aa + a->i[row]);
2220: if (idx) {
2221: itmp = a->j + a->i[row];
2222: if (*nz) *idx = itmp;
2223: else *idx = NULL;
2224: }
2225: MatSeqAIJRestoreArrayRead(A,&aa);
2226: return 0;
2227: }
2229: PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
2230: {
2231: if (nz) *nz = 0;
2232: if (idx) *idx = NULL;
2233: if (v) *v = NULL;
2234: return 0;
2235: }
2237: PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm)
2238: {
2239: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2240: const MatScalar *v;
2241: PetscReal sum = 0.0;
2242: PetscInt i,j;
2244: MatSeqAIJGetArrayRead(A,&v);
2245: if (type == NORM_FROBENIUS) {
2246: #if defined(PETSC_USE_REAL___FP16)
2247: PetscBLASInt one = 1,nz = a->nz;
2248: PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one));
2249: #else
2250: for (i=0; i<a->nz; i++) {
2251: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
2252: }
2253: *nrm = PetscSqrtReal(sum);
2254: #endif
2255: PetscLogFlops(2.0*a->nz);
2256: } else if (type == NORM_1) {
2257: PetscReal *tmp;
2258: PetscInt *jj = a->j;
2259: PetscCalloc1(A->cmap->n+1,&tmp);
2260: *nrm = 0.0;
2261: for (j=0; j<a->nz; j++) {
2262: tmp[*jj++] += PetscAbsScalar(*v); v++;
2263: }
2264: for (j=0; j<A->cmap->n; j++) {
2265: if (tmp[j] > *nrm) *nrm = tmp[j];
2266: }
2267: PetscFree(tmp);
2268: PetscLogFlops(PetscMax(a->nz-1,0));
2269: } else if (type == NORM_INFINITY) {
2270: *nrm = 0.0;
2271: for (j=0; j<A->rmap->n; j++) {
2272: const PetscScalar *v2 = v + a->i[j];
2273: sum = 0.0;
2274: for (i=0; i<a->i[j+1]-a->i[j]; i++) {
2275: sum += PetscAbsScalar(*v2); v2++;
2276: }
2277: if (sum > *nrm) *nrm = sum;
2278: }
2279: PetscLogFlops(PetscMax(a->nz-1,0));
2280: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm");
2281: MatSeqAIJRestoreArrayRead(A,&v);
2282: return 0;
2283: }
2285: /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */
2286: PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B)
2287: {
2288: PetscInt i,j,anzj;
2289: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b;
2290: PetscInt an=A->cmap->N,am=A->rmap->N;
2291: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
2293: /* Allocate space for symbolic transpose info and work array */
2294: PetscCalloc1(an+1,&ati);
2295: PetscMalloc1(ai[am],&atj);
2296: PetscMalloc1(an,&atfill);
2298: /* Walk through aj and count ## of non-zeros in each row of A^T. */
2299: /* Note: offset by 1 for fast conversion into csr format. */
2300: for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1;
2301: /* Form ati for csr format of A^T. */
2302: for (i=0;i<an;i++) ati[i+1] += ati[i];
2304: /* Copy ati into atfill so we have locations of the next free space in atj */
2305: PetscArraycpy(atfill,ati,an);
2307: /* Walk through A row-wise and mark nonzero entries of A^T. */
2308: for (i=0;i<am;i++) {
2309: anzj = ai[i+1] - ai[i];
2310: for (j=0;j<anzj;j++) {
2311: atj[atfill[*aj]] = i;
2312: atfill[*aj++] += 1;
2313: }
2314: }
2316: /* Clean up temporary space and complete requests. */
2317: PetscFree(atfill);
2318: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);
2319: MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
2320: MatSetType(*B,((PetscObject)A)->type_name);
2322: b = (Mat_SeqAIJ*)((*B)->data);
2323: b->free_a = PETSC_FALSE;
2324: b->free_ij = PETSC_TRUE;
2325: b->nonew = 0;
2326: return 0;
2327: }
2329: PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2330: {
2331: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data;
2332: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2333: const MatScalar *va,*vb;
2334: PetscInt ma,na,mb,nb, i;
2336: MatGetSize(A,&ma,&na);
2337: MatGetSize(B,&mb,&nb);
2338: if (ma!=nb || na!=mb) {
2339: *f = PETSC_FALSE;
2340: return 0;
2341: }
2342: MatSeqAIJGetArrayRead(A,&va);
2343: MatSeqAIJGetArrayRead(B,&vb);
2344: aii = aij->i; bii = bij->i;
2345: adx = aij->j; bdx = bij->j;
2346: PetscMalloc1(ma,&aptr);
2347: PetscMalloc1(mb,&bptr);
2348: for (i=0; i<ma; i++) aptr[i] = aii[i];
2349: for (i=0; i<mb; i++) bptr[i] = bii[i];
2351: *f = PETSC_TRUE;
2352: for (i=0; i<ma; i++) {
2353: while (aptr[i]<aii[i+1]) {
2354: PetscInt idc,idr;
2355: PetscScalar vc,vr;
2356: /* column/row index/value */
2357: idc = adx[aptr[i]];
2358: idr = bdx[bptr[idc]];
2359: vc = va[aptr[i]];
2360: vr = vb[bptr[idc]];
2361: if (i!=idr || PetscAbsScalar(vc-vr) > tol) {
2362: *f = PETSC_FALSE;
2363: goto done;
2364: } else {
2365: aptr[i]++;
2366: if (B || i!=idc) bptr[idc]++;
2367: }
2368: }
2369: }
2370: done:
2371: PetscFree(aptr);
2372: PetscFree(bptr);
2373: MatSeqAIJRestoreArrayRead(A,&va);
2374: MatSeqAIJRestoreArrayRead(B,&vb);
2375: return 0;
2376: }
2378: PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2379: {
2380: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data;
2381: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2382: MatScalar *va,*vb;
2383: PetscInt ma,na,mb,nb, i;
2385: MatGetSize(A,&ma,&na);
2386: MatGetSize(B,&mb,&nb);
2387: if (ma!=nb || na!=mb) {
2388: *f = PETSC_FALSE;
2389: return 0;
2390: }
2391: aii = aij->i; bii = bij->i;
2392: adx = aij->j; bdx = bij->j;
2393: va = aij->a; vb = bij->a;
2394: PetscMalloc1(ma,&aptr);
2395: PetscMalloc1(mb,&bptr);
2396: for (i=0; i<ma; i++) aptr[i] = aii[i];
2397: for (i=0; i<mb; i++) bptr[i] = bii[i];
2399: *f = PETSC_TRUE;
2400: for (i=0; i<ma; i++) {
2401: while (aptr[i]<aii[i+1]) {
2402: PetscInt idc,idr;
2403: PetscScalar vc,vr;
2404: /* column/row index/value */
2405: idc = adx[aptr[i]];
2406: idr = bdx[bptr[idc]];
2407: vc = va[aptr[i]];
2408: vr = vb[bptr[idc]];
2409: if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) {
2410: *f = PETSC_FALSE;
2411: goto done;
2412: } else {
2413: aptr[i]++;
2414: if (B || i!=idc) bptr[idc]++;
2415: }
2416: }
2417: }
2418: done:
2419: PetscFree(aptr);
2420: PetscFree(bptr);
2421: return 0;
2422: }
2424: PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2425: {
2426: MatIsTranspose_SeqAIJ(A,A,tol,f);
2427: return 0;
2428: }
2430: PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2431: {
2432: MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);
2433: return 0;
2434: }
2436: PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr)
2437: {
2438: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2439: const PetscScalar *l,*r;
2440: PetscScalar x;
2441: MatScalar *v;
2442: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz;
2443: const PetscInt *jj;
2445: if (ll) {
2446: /* The local size is used so that VecMPI can be passed to this routine
2447: by MatDiagonalScale_MPIAIJ */
2448: VecGetLocalSize(ll,&m);
2450: VecGetArrayRead(ll,&l);
2451: MatSeqAIJGetArray(A,&v);
2452: for (i=0; i<m; i++) {
2453: x = l[i];
2454: M = a->i[i+1] - a->i[i];
2455: for (j=0; j<M; j++) (*v++) *= x;
2456: }
2457: VecRestoreArrayRead(ll,&l);
2458: PetscLogFlops(nz);
2459: MatSeqAIJRestoreArray(A,&v);
2460: }
2461: if (rr) {
2462: VecGetLocalSize(rr,&n);
2464: VecGetArrayRead(rr,&r);
2465: MatSeqAIJGetArray(A,&v);
2466: jj = a->j;
2467: for (i=0; i<nz; i++) (*v++) *= r[*jj++];
2468: MatSeqAIJRestoreArray(A,&v);
2469: VecRestoreArrayRead(rr,&r);
2470: PetscLogFlops(nz);
2471: }
2472: MatSeqAIJInvalidateDiagonal(A);
2473: return 0;
2474: }
2476: PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B)
2477: {
2478: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c;
2479: PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens;
2480: PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi;
2481: const PetscInt *irow,*icol;
2482: const PetscScalar *aa;
2483: PetscInt nrows,ncols;
2484: PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen;
2485: MatScalar *a_new,*mat_a;
2486: Mat C;
2487: PetscBool stride;
2489: ISGetIndices(isrow,&irow);
2490: ISGetLocalSize(isrow,&nrows);
2491: ISGetLocalSize(iscol,&ncols);
2493: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
2494: if (stride) {
2495: ISStrideGetInfo(iscol,&first,&step);
2496: } else {
2497: first = 0;
2498: step = 0;
2499: }
2500: if (stride && step == 1) {
2501: /* special case of contiguous rows */
2502: PetscMalloc2(nrows,&lens,nrows,&starts);
2503: /* loop over new rows determining lens and starting points */
2504: for (i=0; i<nrows; i++) {
2505: kstart = ai[irow[i]];
2506: kend = kstart + ailen[irow[i]];
2507: starts[i] = kstart;
2508: for (k=kstart; k<kend; k++) {
2509: if (aj[k] >= first) {
2510: starts[i] = k;
2511: break;
2512: }
2513: }
2514: sum = 0;
2515: while (k < kend) {
2516: if (aj[k++] >= first+ncols) break;
2517: sum++;
2518: }
2519: lens[i] = sum;
2520: }
2521: /* create submatrix */
2522: if (scall == MAT_REUSE_MATRIX) {
2523: PetscInt n_cols,n_rows;
2524: MatGetSize(*B,&n_rows,&n_cols);
2526: MatZeroEntries(*B);
2527: C = *B;
2528: } else {
2529: PetscInt rbs,cbs;
2530: MatCreate(PetscObjectComm((PetscObject)A),&C);
2531: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2532: ISGetBlockSize(isrow,&rbs);
2533: ISGetBlockSize(iscol,&cbs);
2534: MatSetBlockSizes(C,rbs,cbs);
2535: MatSetType(C,((PetscObject)A)->type_name);
2536: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2537: }
2538: c = (Mat_SeqAIJ*)C->data;
2540: /* loop over rows inserting into submatrix */
2541: a_new = c->a;
2542: j_new = c->j;
2543: i_new = c->i;
2544: MatSeqAIJGetArrayRead(A,&aa);
2545: for (i=0; i<nrows; i++) {
2546: ii = starts[i];
2547: lensi = lens[i];
2548: for (k=0; k<lensi; k++) {
2549: *j_new++ = aj[ii+k] - first;
2550: }
2551: PetscArraycpy(a_new,aa + starts[i],lensi);
2552: a_new += lensi;
2553: i_new[i+1] = i_new[i] + lensi;
2554: c->ilen[i] = lensi;
2555: }
2556: MatSeqAIJRestoreArrayRead(A,&aa);
2557: PetscFree2(lens,starts);
2558: } else {
2559: ISGetIndices(iscol,&icol);
2560: PetscCalloc1(oldcols,&smap);
2561: PetscMalloc1(1+nrows,&lens);
2562: for (i=0; i<ncols; i++) {
2564: smap[icol[i]] = i+1;
2565: }
2567: /* determine lens of each row */
2568: for (i=0; i<nrows; i++) {
2569: kstart = ai[irow[i]];
2570: kend = kstart + a->ilen[irow[i]];
2571: lens[i] = 0;
2572: for (k=kstart; k<kend; k++) {
2573: if (smap[aj[k]]) {
2574: lens[i]++;
2575: }
2576: }
2577: }
2578: /* Create and fill new matrix */
2579: if (scall == MAT_REUSE_MATRIX) {
2580: PetscBool equal;
2582: c = (Mat_SeqAIJ*)((*B)->data);
2584: PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal);
2586: PetscArrayzero(c->ilen,(*B)->rmap->n);
2587: C = *B;
2588: } else {
2589: PetscInt rbs,cbs;
2590: MatCreate(PetscObjectComm((PetscObject)A),&C);
2591: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2592: ISGetBlockSize(isrow,&rbs);
2593: ISGetBlockSize(iscol,&cbs);
2594: MatSetBlockSizes(C,rbs,cbs);
2595: MatSetType(C,((PetscObject)A)->type_name);
2596: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2597: }
2598: MatSeqAIJGetArrayRead(A,&aa);
2599: c = (Mat_SeqAIJ*)(C->data);
2600: for (i=0; i<nrows; i++) {
2601: row = irow[i];
2602: kstart = ai[row];
2603: kend = kstart + a->ilen[row];
2604: mat_i = c->i[i];
2605: mat_j = c->j + mat_i;
2606: mat_a = c->a + mat_i;
2607: mat_ilen = c->ilen + i;
2608: for (k=kstart; k<kend; k++) {
2609: if ((tcol=smap[a->j[k]])) {
2610: *mat_j++ = tcol - 1;
2611: *mat_a++ = aa[k];
2612: (*mat_ilen)++;
2614: }
2615: }
2616: }
2617: MatSeqAIJRestoreArrayRead(A,&aa);
2618: /* Free work space */
2619: ISRestoreIndices(iscol,&icol);
2620: PetscFree(smap);
2621: PetscFree(lens);
2622: /* sort */
2623: for (i = 0; i < nrows; i++) {
2624: PetscInt ilen;
2626: mat_i = c->i[i];
2627: mat_j = c->j + mat_i;
2628: mat_a = c->a + mat_i;
2629: ilen = c->ilen[i];
2630: PetscSortIntWithScalarArray(ilen,mat_j,mat_a);
2631: }
2632: }
2633: #if defined(PETSC_HAVE_DEVICE)
2634: MatBindToCPU(C,A->boundtocpu);
2635: #endif
2636: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2637: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2639: ISRestoreIndices(isrow,&irow);
2640: *B = C;
2641: return 0;
2642: }
2644: PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat)
2645: {
2646: Mat B;
2648: if (scall == MAT_INITIAL_MATRIX) {
2649: MatCreate(subComm,&B);
2650: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);
2651: MatSetBlockSizesFromMats(B,mat,mat);
2652: MatSetType(B,MATSEQAIJ);
2653: MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);
2654: *subMat = B;
2655: } else {
2656: MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);
2657: }
2658: return 0;
2659: }
2661: PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2662: {
2663: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2664: Mat outA;
2665: PetscBool row_identity,col_identity;
2669: ISIdentity(row,&row_identity);
2670: ISIdentity(col,&col_identity);
2672: outA = inA;
2673: outA->factortype = MAT_FACTOR_LU;
2674: PetscFree(inA->solvertype);
2675: PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);
2677: PetscObjectReference((PetscObject)row);
2678: ISDestroy(&a->row);
2680: a->row = row;
2682: PetscObjectReference((PetscObject)col);
2683: ISDestroy(&a->col);
2685: a->col = col;
2687: /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */
2688: ISDestroy(&a->icol);
2689: ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2690: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
2692: if (!a->solve_work) { /* this matrix may have been factored before */
2693: PetscMalloc1(inA->rmap->n+1,&a->solve_work);
2694: PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));
2695: }
2697: MatMarkDiagonal_SeqAIJ(inA);
2698: if (row_identity && col_identity) {
2699: MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);
2700: } else {
2701: MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);
2702: }
2703: return 0;
2704: }
2706: PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha)
2707: {
2708: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2709: PetscScalar *v;
2710: PetscBLASInt one = 1,bnz;
2712: MatSeqAIJGetArray(inA,&v);
2713: PetscBLASIntCast(a->nz,&bnz);
2714: PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one));
2715: PetscLogFlops(a->nz);
2716: MatSeqAIJRestoreArray(inA,&v);
2717: MatSeqAIJInvalidateDiagonal(inA);
2718: return 0;
2719: }
2721: PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj)
2722: {
2723: PetscInt i;
2725: if (!submatj->id) { /* delete data that are linked only to submats[id=0] */
2726: PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);
2728: for (i=0; i<submatj->nrqr; ++i) {
2729: PetscFree(submatj->sbuf2[i]);
2730: }
2731: PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);
2733: if (submatj->rbuf1) {
2734: PetscFree(submatj->rbuf1[0]);
2735: PetscFree(submatj->rbuf1);
2736: }
2738: for (i=0; i<submatj->nrqs; ++i) {
2739: PetscFree(submatj->rbuf3[i]);
2740: }
2741: PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);
2742: PetscFree(submatj->pa);
2743: }
2745: #if defined(PETSC_USE_CTABLE)
2746: PetscTableDestroy((PetscTable*)&submatj->rmap);
2747: if (submatj->cmap_loc) PetscFree(submatj->cmap_loc);
2748: PetscFree(submatj->rmap_loc);
2749: #else
2750: PetscFree(submatj->rmap);
2751: #endif
2753: if (!submatj->allcolumns) {
2754: #if defined(PETSC_USE_CTABLE)
2755: PetscTableDestroy((PetscTable*)&submatj->cmap);
2756: #else
2757: PetscFree(submatj->cmap);
2758: #endif
2759: }
2760: PetscFree(submatj->row2proc);
2762: PetscFree(submatj);
2763: return 0;
2764: }
2766: PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C)
2767: {
2768: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
2769: Mat_SubSppt *submatj = c->submatis1;
2771: (*submatj->destroy)(C);
2772: MatDestroySubMatrix_Private(submatj);
2773: return 0;
2774: }
2776: PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[])
2777: {
2778: PetscInt i;
2779: Mat C;
2780: Mat_SeqAIJ *c;
2781: Mat_SubSppt *submatj;
2783: for (i=0; i<n; i++) {
2784: C = (*mat)[i];
2785: c = (Mat_SeqAIJ*)C->data;
2786: submatj = c->submatis1;
2787: if (submatj) {
2788: if (--((PetscObject)C)->refct <= 0) {
2789: (*submatj->destroy)(C);
2790: MatDestroySubMatrix_Private(submatj);
2791: PetscFree(C->defaultvectype);
2792: PetscLayoutDestroy(&C->rmap);
2793: PetscLayoutDestroy(&C->cmap);
2794: PetscHeaderDestroy(&C);
2795: }
2796: } else {
2797: MatDestroy(&C);
2798: }
2799: }
2801: /* Destroy Dummy submatrices created for reuse */
2802: MatDestroySubMatrices_Dummy(n,mat);
2804: PetscFree(*mat);
2805: return 0;
2806: }
2808: PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2809: {
2810: PetscInt i;
2812: if (scall == MAT_INITIAL_MATRIX) {
2813: PetscCalloc1(n+1,B);
2814: }
2816: for (i=0; i<n; i++) {
2817: MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2818: }
2819: return 0;
2820: }
2822: PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov)
2823: {
2824: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2825: PetscInt row,i,j,k,l,m,n,*nidx,isz,val;
2826: const PetscInt *idx;
2827: PetscInt start,end,*ai,*aj;
2828: PetscBT table;
2830: m = A->rmap->n;
2831: ai = a->i;
2832: aj = a->j;
2836: PetscMalloc1(m+1,&nidx);
2837: PetscBTCreate(m,&table);
2839: for (i=0; i<is_max; i++) {
2840: /* Initialize the two local arrays */
2841: isz = 0;
2842: PetscBTMemzero(m,table);
2844: /* Extract the indices, assume there can be duplicate entries */
2845: ISGetIndices(is[i],&idx);
2846: ISGetLocalSize(is[i],&n);
2848: /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */
2849: for (j=0; j<n; ++j) {
2850: if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j];
2851: }
2852: ISRestoreIndices(is[i],&idx);
2853: ISDestroy(&is[i]);
2855: k = 0;
2856: for (j=0; j<ov; j++) { /* for each overlap */
2857: n = isz;
2858: for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */
2859: row = nidx[k];
2860: start = ai[row];
2861: end = ai[row+1];
2862: for (l = start; l<end; l++) {
2863: val = aj[l];
2864: if (!PetscBTLookupSet(table,val)) nidx[isz++] = val;
2865: }
2866: }
2867: }
2868: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));
2869: }
2870: PetscBTDestroy(&table);
2871: PetscFree(nidx);
2872: return 0;
2873: }
2875: /* -------------------------------------------------------------- */
2876: PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B)
2877: {
2878: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2879: PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n;
2880: const PetscInt *row,*col;
2881: PetscInt *cnew,j,*lens;
2882: IS icolp,irowp;
2883: PetscInt *cwork = NULL;
2884: PetscScalar *vwork = NULL;
2886: ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);
2887: ISGetIndices(irowp,&row);
2888: ISInvertPermutation(colp,PETSC_DECIDE,&icolp);
2889: ISGetIndices(icolp,&col);
2891: /* determine lengths of permuted rows */
2892: PetscMalloc1(m+1,&lens);
2893: for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i];
2894: MatCreate(PetscObjectComm((PetscObject)A),B);
2895: MatSetSizes(*B,m,n,m,n);
2896: MatSetBlockSizesFromMats(*B,A,A);
2897: MatSetType(*B,((PetscObject)A)->type_name);
2898: MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);
2899: PetscFree(lens);
2901: PetscMalloc1(n,&cnew);
2902: for (i=0; i<m; i++) {
2903: MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2904: for (j=0; j<nz; j++) cnew[j] = col[cwork[j]];
2905: MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);
2906: MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2907: }
2908: PetscFree(cnew);
2910: (*B)->assembled = PETSC_FALSE;
2912: #if defined(PETSC_HAVE_DEVICE)
2913: MatBindToCPU(*B,A->boundtocpu);
2914: #endif
2915: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
2916: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
2917: ISRestoreIndices(irowp,&row);
2918: ISRestoreIndices(icolp,&col);
2919: ISDestroy(&irowp);
2920: ISDestroy(&icolp);
2921: if (rowp == colp) {
2922: MatPropagateSymmetryOptions(A,*B);
2923: }
2924: return 0;
2925: }
2927: PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str)
2928: {
2929: /* If the two matrices have the same copy implementation, use fast copy. */
2930: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2931: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2932: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
2933: const PetscScalar *aa;
2935: MatSeqAIJGetArrayRead(A,&aa);
2937: PetscArraycpy(b->a,aa,a->i[A->rmap->n]);
2938: PetscObjectStateIncrease((PetscObject)B);
2939: MatSeqAIJRestoreArrayRead(A,&aa);
2940: } else {
2941: MatCopy_Basic(A,B,str);
2942: }
2943: return 0;
2944: }
2946: PetscErrorCode MatSetUp_SeqAIJ(Mat A)
2947: {
2948: MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL);
2949: return 0;
2950: }
2952: PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[])
2953: {
2954: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2956: *array = a->a;
2957: return 0;
2958: }
2960: PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[])
2961: {
2962: *array = NULL;
2963: return 0;
2964: }
2966: /*
2967: Computes the number of nonzeros per row needed for preallocation when X and Y
2968: have different nonzero structure.
2969: */
2970: PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz)
2971: {
2972: PetscInt i,j,k,nzx,nzy;
2974: /* Set the number of nonzeros in the new matrix */
2975: for (i=0; i<m; i++) {
2976: const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i];
2977: nzx = xi[i+1] - xi[i];
2978: nzy = yi[i+1] - yi[i];
2979: nnz[i] = 0;
2980: for (j=0,k=0; j<nzx; j++) { /* Point in X */
2981: for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */
2982: if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */
2983: nnz[i]++;
2984: }
2985: for (; k<nzy; k++) nnz[i]++;
2986: }
2987: return 0;
2988: }
2990: PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz)
2991: {
2992: PetscInt m = Y->rmap->N;
2993: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data;
2994: Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data;
2996: /* Set the number of nonzeros in the new matrix */
2997: MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);
2998: return 0;
2999: }
3001: PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
3002: {
3003: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data;
3005: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
3006: PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE;
3007: if (e) {
3008: PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e);
3009: if (e) {
3010: PetscArraycmp(x->j,y->j,y->nz,&e);
3011: if (e) str = SAME_NONZERO_PATTERN;
3012: }
3013: }
3015: }
3016: if (str == SAME_NONZERO_PATTERN) {
3017: const PetscScalar *xa;
3018: PetscScalar *ya,alpha = a;
3019: PetscBLASInt one = 1,bnz;
3021: PetscBLASIntCast(x->nz,&bnz);
3022: MatSeqAIJGetArray(Y,&ya);
3023: MatSeqAIJGetArrayRead(X,&xa);
3024: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one));
3025: MatSeqAIJRestoreArrayRead(X,&xa);
3026: MatSeqAIJRestoreArray(Y,&ya);
3027: PetscLogFlops(2.0*bnz);
3028: MatSeqAIJInvalidateDiagonal(Y);
3029: PetscObjectStateIncrease((PetscObject)Y);
3030: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
3031: MatAXPY_Basic(Y,a,X,str);
3032: } else {
3033: Mat B;
3034: PetscInt *nnz;
3035: PetscMalloc1(Y->rmap->N,&nnz);
3036: MatCreate(PetscObjectComm((PetscObject)Y),&B);
3037: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
3038: MatSetLayouts(B,Y->rmap,Y->cmap);
3039: MatSetType(B,((PetscObject)Y)->type_name);
3040: MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);
3041: MatSeqAIJSetPreallocation(B,0,nnz);
3042: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
3043: MatHeaderMerge(Y,&B);
3044: PetscFree(nnz);
3045: }
3046: return 0;
3047: }
3049: PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat)
3050: {
3051: #if defined(PETSC_USE_COMPLEX)
3052: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3053: PetscInt i,nz;
3054: PetscScalar *a;
3056: nz = aij->nz;
3057: MatSeqAIJGetArray(mat,&a);
3058: for (i=0; i<nz; i++) a[i] = PetscConj(a[i]);
3059: MatSeqAIJRestoreArray(mat,&a);
3060: #else
3061: #endif
3062: return 0;
3063: }
3065: PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3066: {
3067: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3068: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3069: PetscReal atmp;
3070: PetscScalar *x;
3071: const MatScalar *aa,*av;
3074: MatSeqAIJGetArrayRead(A,&av);
3075: aa = av;
3076: ai = a->i;
3077: aj = a->j;
3079: VecSet(v,0.0);
3080: VecGetArrayWrite(v,&x);
3081: VecGetLocalSize(v,&n);
3083: for (i=0; i<m; i++) {
3084: ncols = ai[1] - ai[0]; ai++;
3085: for (j=0; j<ncols; j++) {
3086: atmp = PetscAbsScalar(*aa);
3087: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
3088: aa++; aj++;
3089: }
3090: }
3091: VecRestoreArrayWrite(v,&x);
3092: MatSeqAIJRestoreArrayRead(A,&av);
3093: return 0;
3094: }
3096: PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3097: {
3098: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3099: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3100: PetscScalar *x;
3101: const MatScalar *aa,*av;
3104: MatSeqAIJGetArrayRead(A,&av);
3105: aa = av;
3106: ai = a->i;
3107: aj = a->j;
3109: VecSet(v,0.0);
3110: VecGetArrayWrite(v,&x);
3111: VecGetLocalSize(v,&n);
3113: for (i=0; i<m; i++) {
3114: ncols = ai[1] - ai[0]; ai++;
3115: if (ncols == A->cmap->n) { /* row is dense */
3116: x[i] = *aa; if (idx) idx[i] = 0;
3117: } else { /* row is sparse so already KNOW maximum is 0.0 or higher */
3118: x[i] = 0.0;
3119: if (idx) {
3120: for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */
3121: if (aj[j] > j) {
3122: idx[i] = j;
3123: break;
3124: }
3125: }
3126: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3127: if (j==ncols && j < A->cmap->n) idx[i] = j;
3128: }
3129: }
3130: for (j=0; j<ncols; j++) {
3131: if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3132: aa++; aj++;
3133: }
3134: }
3135: VecRestoreArrayWrite(v,&x);
3136: MatSeqAIJRestoreArrayRead(A,&av);
3137: return 0;
3138: }
3140: PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3141: {
3142: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3143: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3144: PetscScalar *x;
3145: const MatScalar *aa,*av;
3147: MatSeqAIJGetArrayRead(A,&av);
3148: aa = av;
3149: ai = a->i;
3150: aj = a->j;
3152: VecSet(v,0.0);
3153: VecGetArrayWrite(v,&x);
3154: VecGetLocalSize(v,&n);
3156: for (i=0; i<m; i++) {
3157: ncols = ai[1] - ai[0]; ai++;
3158: if (ncols == A->cmap->n) { /* row is dense */
3159: x[i] = *aa; if (idx) idx[i] = 0;
3160: } else { /* row is sparse so already KNOW minimum is 0.0 or higher */
3161: x[i] = 0.0;
3162: if (idx) { /* find first implicit 0.0 in the row */
3163: for (j=0; j<ncols; j++) {
3164: if (aj[j] > j) {
3165: idx[i] = j;
3166: break;
3167: }
3168: }
3169: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3170: if (j==ncols && j < A->cmap->n) idx[i] = j;
3171: }
3172: }
3173: for (j=0; j<ncols; j++) {
3174: if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3175: aa++; aj++;
3176: }
3177: }
3178: VecRestoreArrayWrite(v,&x);
3179: MatSeqAIJRestoreArrayRead(A,&av);
3180: return 0;
3181: }
3183: PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3184: {
3185: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3186: PetscInt i,j,m = A->rmap->n,ncols,n;
3187: const PetscInt *ai,*aj;
3188: PetscScalar *x;
3189: const MatScalar *aa,*av;
3192: MatSeqAIJGetArrayRead(A,&av);
3193: aa = av;
3194: ai = a->i;
3195: aj = a->j;
3197: VecSet(v,0.0);
3198: VecGetArrayWrite(v,&x);
3199: VecGetLocalSize(v,&n);
3201: for (i=0; i<m; i++) {
3202: ncols = ai[1] - ai[0]; ai++;
3203: if (ncols == A->cmap->n) { /* row is dense */
3204: x[i] = *aa; if (idx) idx[i] = 0;
3205: } else { /* row is sparse so already KNOW minimum is 0.0 or lower */
3206: x[i] = 0.0;
3207: if (idx) { /* find first implicit 0.0 in the row */
3208: for (j=0; j<ncols; j++) {
3209: if (aj[j] > j) {
3210: idx[i] = j;
3211: break;
3212: }
3213: }
3214: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3215: if (j==ncols && j < A->cmap->n) idx[i] = j;
3216: }
3217: }
3218: for (j=0; j<ncols; j++) {
3219: if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3220: aa++; aj++;
3221: }
3222: }
3223: VecRestoreArrayWrite(v,&x);
3224: MatSeqAIJRestoreArrayRead(A,&av);
3225: return 0;
3226: }
3228: PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values)
3229: {
3230: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
3231: PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j;
3232: MatScalar *diag,work[25],*v_work;
3233: const PetscReal shift = 0.0;
3234: PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
3236: allowzeropivot = PetscNot(A->erroriffailure);
3237: if (a->ibdiagvalid) {
3238: if (values) *values = a->ibdiag;
3239: return 0;
3240: }
3241: MatMarkDiagonal_SeqAIJ(A);
3242: if (!a->ibdiag) {
3243: PetscMalloc1(bs2*mbs,&a->ibdiag);
3244: PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));
3245: }
3246: diag = a->ibdiag;
3247: if (values) *values = a->ibdiag;
3248: /* factor and invert each block */
3249: switch (bs) {
3250: case 1:
3251: for (i=0; i<mbs; i++) {
3252: MatGetValues(A,1,&i,1,&i,diag+i);
3253: if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) {
3254: if (allowzeropivot) {
3255: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3256: A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]);
3257: A->factorerror_zeropivot_row = i;
3258: PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);
3259: } 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);
3260: }
3261: diag[i] = (PetscScalar)1.0 / (diag[i] + shift);
3262: }
3263: break;
3264: case 2:
3265: for (i=0; i<mbs; i++) {
3266: ij[0] = 2*i; ij[1] = 2*i + 1;
3267: MatGetValues(A,2,ij,2,ij,diag);
3268: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
3269: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3270: PetscKernel_A_gets_transpose_A_2(diag);
3271: diag += 4;
3272: }
3273: break;
3274: case 3:
3275: for (i=0; i<mbs; i++) {
3276: ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2;
3277: MatGetValues(A,3,ij,3,ij,diag);
3278: PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
3279: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3280: PetscKernel_A_gets_transpose_A_3(diag);
3281: diag += 9;
3282: }
3283: break;
3284: case 4:
3285: for (i=0; i<mbs; i++) {
3286: ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3;
3287: MatGetValues(A,4,ij,4,ij,diag);
3288: PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
3289: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3290: PetscKernel_A_gets_transpose_A_4(diag);
3291: diag += 16;
3292: }
3293: break;
3294: case 5:
3295: for (i=0; i<mbs; i++) {
3296: ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4;
3297: MatGetValues(A,5,ij,5,ij,diag);
3298: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
3299: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3300: PetscKernel_A_gets_transpose_A_5(diag);
3301: diag += 25;
3302: }
3303: break;
3304: case 6:
3305: for (i=0; i<mbs; i++) {
3306: 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;
3307: MatGetValues(A,6,ij,6,ij,diag);
3308: PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
3309: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3310: PetscKernel_A_gets_transpose_A_6(diag);
3311: diag += 36;
3312: }
3313: break;
3314: case 7:
3315: for (i=0; i<mbs; i++) {
3316: 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;
3317: MatGetValues(A,7,ij,7,ij,diag);
3318: PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
3319: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3320: PetscKernel_A_gets_transpose_A_7(diag);
3321: diag += 49;
3322: }
3323: break;
3324: default:
3325: PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);
3326: for (i=0; i<mbs; i++) {
3327: for (j=0; j<bs; j++) {
3328: IJ[j] = bs*i + j;
3329: }
3330: MatGetValues(A,bs,IJ,bs,IJ,diag);
3331: PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
3332: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3333: PetscKernel_A_gets_transpose_A_N(diag,bs);
3334: diag += bs2;
3335: }
3336: PetscFree3(v_work,v_pivots,IJ);
3337: }
3338: a->ibdiagvalid = PETSC_TRUE;
3339: return 0;
3340: }
3342: static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx)
3343: {
3344: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3345: PetscScalar a,*aa;
3346: PetscInt m,n,i,j,col;
3348: if (!x->assembled) {
3349: MatGetSize(x,&m,&n);
3350: for (i=0; i<m; i++) {
3351: for (j=0; j<aij->imax[i]; j++) {
3352: PetscRandomGetValue(rctx,&a);
3353: col = (PetscInt)(n*PetscRealPart(a));
3354: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3355: }
3356: }
3357: } else {
3358: MatSeqAIJGetArrayWrite(x,&aa);
3359: for (i=0; i<aij->nz; i++) PetscRandomGetValue(rctx,aa+i);
3360: MatSeqAIJRestoreArrayWrite(x,&aa);
3361: }
3362: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3363: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3364: return 0;
3365: }
3367: /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */
3368: PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx)
3369: {
3370: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3371: PetscScalar a;
3372: PetscInt m,n,i,j,col,nskip;
3374: nskip = high - low;
3375: MatGetSize(x,&m,&n);
3376: n -= nskip; /* shrink number of columns where nonzeros can be set */
3377: for (i=0; i<m; i++) {
3378: for (j=0; j<aij->imax[i]; j++) {
3379: PetscRandomGetValue(rctx,&a);
3380: col = (PetscInt)(n*PetscRealPart(a));
3381: if (col >= low) col += nskip; /* shift col rightward to skip the hole */
3382: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3383: }
3384: }
3385: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3386: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3387: return 0;
3388: }
3390: /* -------------------------------------------------------------------*/
3391: static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ,
3392: MatGetRow_SeqAIJ,
3393: MatRestoreRow_SeqAIJ,
3394: MatMult_SeqAIJ,
3395: /* 4*/ MatMultAdd_SeqAIJ,
3396: MatMultTranspose_SeqAIJ,
3397: MatMultTransposeAdd_SeqAIJ,
3398: NULL,
3399: NULL,
3400: NULL,
3401: /* 10*/ NULL,
3402: MatLUFactor_SeqAIJ,
3403: NULL,
3404: MatSOR_SeqAIJ,
3405: MatTranspose_SeqAIJ,
3406: /*1 5*/ MatGetInfo_SeqAIJ,
3407: MatEqual_SeqAIJ,
3408: MatGetDiagonal_SeqAIJ,
3409: MatDiagonalScale_SeqAIJ,
3410: MatNorm_SeqAIJ,
3411: /* 20*/ NULL,
3412: MatAssemblyEnd_SeqAIJ,
3413: MatSetOption_SeqAIJ,
3414: MatZeroEntries_SeqAIJ,
3415: /* 24*/ MatZeroRows_SeqAIJ,
3416: NULL,
3417: NULL,
3418: NULL,
3419: NULL,
3420: /* 29*/ MatSetUp_SeqAIJ,
3421: NULL,
3422: NULL,
3423: NULL,
3424: NULL,
3425: /* 34*/ MatDuplicate_SeqAIJ,
3426: NULL,
3427: NULL,
3428: MatILUFactor_SeqAIJ,
3429: NULL,
3430: /* 39*/ MatAXPY_SeqAIJ,
3431: MatCreateSubMatrices_SeqAIJ,
3432: MatIncreaseOverlap_SeqAIJ,
3433: MatGetValues_SeqAIJ,
3434: MatCopy_SeqAIJ,
3435: /* 44*/ MatGetRowMax_SeqAIJ,
3436: MatScale_SeqAIJ,
3437: MatShift_SeqAIJ,
3438: MatDiagonalSet_SeqAIJ,
3439: MatZeroRowsColumns_SeqAIJ,
3440: /* 49*/ MatSetRandom_SeqAIJ,
3441: MatGetRowIJ_SeqAIJ,
3442: MatRestoreRowIJ_SeqAIJ,
3443: MatGetColumnIJ_SeqAIJ,
3444: MatRestoreColumnIJ_SeqAIJ,
3445: /* 54*/ MatFDColoringCreate_SeqXAIJ,
3446: NULL,
3447: NULL,
3448: MatPermute_SeqAIJ,
3449: NULL,
3450: /* 59*/ NULL,
3451: MatDestroy_SeqAIJ,
3452: MatView_SeqAIJ,
3453: NULL,
3454: NULL,
3455: /* 64*/ NULL,
3456: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ,
3457: NULL,
3458: NULL,
3459: NULL,
3460: /* 69*/ MatGetRowMaxAbs_SeqAIJ,
3461: MatGetRowMinAbs_SeqAIJ,
3462: NULL,
3463: NULL,
3464: NULL,
3465: /* 74*/ NULL,
3466: MatFDColoringApply_AIJ,
3467: NULL,
3468: NULL,
3469: NULL,
3470: /* 79*/ MatFindZeroDiagonals_SeqAIJ,
3471: NULL,
3472: NULL,
3473: NULL,
3474: MatLoad_SeqAIJ,
3475: /* 84*/ MatIsSymmetric_SeqAIJ,
3476: MatIsHermitian_SeqAIJ,
3477: NULL,
3478: NULL,
3479: NULL,
3480: /* 89*/ NULL,
3481: NULL,
3482: MatMatMultNumeric_SeqAIJ_SeqAIJ,
3483: NULL,
3484: NULL,
3485: /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy,
3486: NULL,
3487: NULL,
3488: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ,
3489: NULL,
3490: /* 99*/ MatProductSetFromOptions_SeqAIJ,
3491: NULL,
3492: NULL,
3493: MatConjugate_SeqAIJ,
3494: NULL,
3495: /*104*/ MatSetValuesRow_SeqAIJ,
3496: MatRealPart_SeqAIJ,
3497: MatImaginaryPart_SeqAIJ,
3498: NULL,
3499: NULL,
3500: /*109*/ MatMatSolve_SeqAIJ,
3501: NULL,
3502: MatGetRowMin_SeqAIJ,
3503: NULL,
3504: MatMissingDiagonal_SeqAIJ,
3505: /*114*/ NULL,
3506: NULL,
3507: NULL,
3508: NULL,
3509: NULL,
3510: /*119*/ NULL,
3511: NULL,
3512: NULL,
3513: NULL,
3514: MatGetMultiProcBlock_SeqAIJ,
3515: /*124*/ MatFindNonzeroRows_SeqAIJ,
3516: MatGetColumnReductions_SeqAIJ,
3517: MatInvertBlockDiagonal_SeqAIJ,
3518: MatInvertVariableBlockDiagonal_SeqAIJ,
3519: NULL,
3520: /*129*/ NULL,
3521: NULL,
3522: NULL,
3523: MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ,
3524: MatTransposeColoringCreate_SeqAIJ,
3525: /*134*/ MatTransColoringApplySpToDen_SeqAIJ,
3526: MatTransColoringApplyDenToSp_SeqAIJ,
3527: NULL,
3528: NULL,
3529: MatRARtNumeric_SeqAIJ_SeqAIJ,
3530: /*139*/NULL,
3531: NULL,
3532: NULL,
3533: MatFDColoringSetUp_SeqXAIJ,
3534: MatFindOffBlockDiagonalEntries_SeqAIJ,
3535: MatCreateMPIMatConcatenateSeqMat_SeqAIJ,
3536: /*145*/MatDestroySubMatrices_SeqAIJ,
3537: NULL,
3538: NULL
3539: };
3541: PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices)
3542: {
3543: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3544: PetscInt i,nz,n;
3546: nz = aij->maxnz;
3547: n = mat->rmap->n;
3548: for (i=0; i<nz; i++) {
3549: aij->j[i] = indices[i];
3550: }
3551: aij->nz = nz;
3552: for (i=0; i<n; i++) {
3553: aij->ilen[i] = aij->imax[i];
3554: }
3555: return 0;
3556: }
3558: /*
3559: * Given a sparse matrix with global column indices, compact it by using a local column space.
3560: * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable()
3561: */
3562: PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping)
3563: {
3564: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3565: PetscTable gid1_lid1;
3566: PetscTablePosition tpos;
3567: PetscInt gid,lid,i,ec,nz = aij->nz;
3568: PetscInt *garray,*jj = aij->j;
3572: /* use a table */
3573: PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);
3574: ec = 0;
3575: for (i=0; i<nz; i++) {
3576: PetscInt data,gid1 = jj[i] + 1;
3577: PetscTableFind(gid1_lid1,gid1,&data);
3578: if (!data) {
3579: /* one based table */
3580: PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);
3581: }
3582: }
3583: /* form array of columns we need */
3584: PetscMalloc1(ec,&garray);
3585: PetscTableGetHeadPosition(gid1_lid1,&tpos);
3586: while (tpos) {
3587: PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);
3588: gid--;
3589: lid--;
3590: garray[lid] = gid;
3591: }
3592: PetscSortInt(ec,garray); /* sort, and rebuild */
3593: PetscTableRemoveAll(gid1_lid1);
3594: for (i=0; i<ec; i++) {
3595: PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);
3596: }
3597: /* compact out the extra columns in B */
3598: for (i=0; i<nz; i++) {
3599: PetscInt gid1 = jj[i] + 1;
3600: PetscTableFind(gid1_lid1,gid1,&lid);
3601: lid--;
3602: jj[i] = lid;
3603: }
3604: PetscLayoutDestroy(&mat->cmap);
3605: PetscTableDestroy(&gid1_lid1);
3606: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap);
3607: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);
3608: ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);
3609: return 0;
3610: }
3612: /*@
3613: MatSeqAIJSetColumnIndices - Set the column indices for all the rows
3614: in the matrix.
3616: Input Parameters:
3617: + mat - the SeqAIJ matrix
3618: - indices - the column indices
3620: Level: advanced
3622: Notes:
3623: This can be called if you have precomputed the nonzero structure of the
3624: matrix and want to provide it to the matrix object to improve the performance
3625: of the MatSetValues() operation.
3627: You MUST have set the correct numbers of nonzeros per row in the call to
3628: MatCreateSeqAIJ(), and the columns indices MUST be sorted.
3630: MUST be called before any calls to MatSetValues();
3632: The indices should start with zero, not one.
3634: @*/
3635: PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices)
3636: {
3639: PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
3640: return 0;
3641: }
3643: /* ----------------------------------------------------------------------------------------*/
3645: PetscErrorCode MatStoreValues_SeqAIJ(Mat mat)
3646: {
3647: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3648: size_t nz = aij->i[mat->rmap->n];
3652: /* allocate space for values if not already there */
3653: if (!aij->saved_values) {
3654: PetscMalloc1(nz+1,&aij->saved_values);
3655: PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
3656: }
3658: /* copy values over */
3659: PetscArraycpy(aij->saved_values,aij->a,nz);
3660: return 0;
3661: }
3663: /*@
3664: MatStoreValues - Stashes a copy of the matrix values; this allows, for
3665: example, reuse of the linear part of a Jacobian, while recomputing the
3666: nonlinear portion.
3668: Collect on Mat
3670: Input Parameters:
3671: . mat - the matrix (currently only AIJ matrices support this option)
3673: Level: advanced
3675: Common Usage, with SNESSolve():
3676: $ Create Jacobian matrix
3677: $ Set linear terms into matrix
3678: $ Apply boundary conditions to matrix, at this time matrix must have
3679: $ final nonzero structure (i.e. setting the nonlinear terms and applying
3680: $ boundary conditions again will not change the nonzero structure
3681: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3682: $ MatStoreValues(mat);
3683: $ Call SNESSetJacobian() with matrix
3684: $ In your Jacobian routine
3685: $ MatRetrieveValues(mat);
3686: $ Set nonlinear terms in matrix
3688: Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself:
3689: $ // build linear portion of Jacobian
3690: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3691: $ MatStoreValues(mat);
3692: $ loop over nonlinear iterations
3693: $ MatRetrieveValues(mat);
3694: $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian
3695: $ // call MatAssemblyBegin/End() on matrix
3696: $ Solve linear system with Jacobian
3697: $ endloop
3699: Notes:
3700: Matrix must already be assemblied before calling this routine
3701: Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before
3702: calling this routine.
3704: When this is called multiple times it overwrites the previous set of stored values
3705: and does not allocated additional space.
3707: .seealso: MatRetrieveValues()
3709: @*/
3710: PetscErrorCode MatStoreValues(Mat mat)
3711: {
3715: PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));
3716: return 0;
3717: }
3719: PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat)
3720: {
3721: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3722: PetscInt nz = aij->i[mat->rmap->n];
3726: /* copy values over */
3727: PetscArraycpy(aij->a,aij->saved_values,nz);
3728: return 0;
3729: }
3731: /*@
3732: MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for
3733: example, reuse of the linear part of a Jacobian, while recomputing the
3734: nonlinear portion.
3736: Collect on Mat
3738: Input Parameters:
3739: . mat - the matrix (currently only AIJ matrices support this option)
3741: Level: advanced
3743: .seealso: MatStoreValues()
3745: @*/
3746: PetscErrorCode MatRetrieveValues(Mat mat)
3747: {
3751: PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));
3752: return 0;
3753: }
3755: /* --------------------------------------------------------------------------------*/
3756: /*@C
3757: MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format
3758: (the default parallel PETSc format). For good matrix assembly performance
3759: the user should preallocate the matrix storage by setting the parameter nz
3760: (or the array nnz). By setting these parameters accurately, performance
3761: during matrix assembly can be increased by more than a factor of 50.
3763: Collective
3765: Input Parameters:
3766: + comm - MPI communicator, set to PETSC_COMM_SELF
3767: . m - number of rows
3768: . n - number of columns
3769: . nz - number of nonzeros per row (same for all rows)
3770: - nnz - array containing the number of nonzeros in the various rows
3771: (possibly different for each row) or NULL
3773: Output Parameter:
3774: . A - the matrix
3776: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3777: MatXXXXSetPreallocation() paradigm instead of this routine directly.
3778: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
3780: Notes:
3781: If nnz is given then nz is ignored
3783: The AIJ format (also called the Yale sparse matrix format or
3784: compressed row storage), is fully compatible with standard Fortran 77
3785: storage. That is, the stored row and column indices can begin at
3786: either one (as in Fortran) or zero. See the users' manual for details.
3788: Specify the preallocated storage with either nz or nnz (not both).
3789: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3790: allocation. For large problems you MUST preallocate memory or you
3791: will get TERRIBLE performance, see the users' manual chapter on matrices.
3793: By default, this format uses inodes (identical nodes) when possible, to
3794: improve numerical efficiency of matrix-vector products and solves. We
3795: search for consecutive rows with the same nonzero structure, thereby
3796: reusing matrix information to achieve increased efficiency.
3798: Options Database Keys:
3799: + -mat_no_inode - Do not use inodes
3800: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3802: Level: intermediate
3804: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays()
3806: @*/
3807: PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3808: {
3809: MatCreate(comm,A);
3810: MatSetSizes(*A,m,n,m,n);
3811: MatSetType(*A,MATSEQAIJ);
3812: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
3813: return 0;
3814: }
3816: /*@C
3817: MatSeqAIJSetPreallocation - For good matrix assembly performance
3818: the user should preallocate the matrix storage by setting the parameter nz
3819: (or the array nnz). By setting these parameters accurately, performance
3820: during matrix assembly can be increased by more than a factor of 50.
3822: Collective
3824: Input Parameters:
3825: + B - The matrix
3826: . nz - number of nonzeros per row (same for all rows)
3827: - nnz - array containing the number of nonzeros in the various rows
3828: (possibly different for each row) or NULL
3830: Notes:
3831: If nnz is given then nz is ignored
3833: The AIJ format (also called the Yale sparse matrix format or
3834: compressed row storage), is fully compatible with standard Fortran 77
3835: storage. That is, the stored row and column indices can begin at
3836: either one (as in Fortran) or zero. See the users' manual for details.
3838: Specify the preallocated storage with either nz or nnz (not both).
3839: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3840: allocation. For large problems you MUST preallocate memory or you
3841: will get TERRIBLE performance, see the users' manual chapter on matrices.
3843: You can call MatGetInfo() to get information on how effective the preallocation was;
3844: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3845: You can also run with the option -info and look for messages with the string
3846: malloc in them to see if additional memory allocation was needed.
3848: Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix
3849: entries or columns indices
3851: By default, this format uses inodes (identical nodes) when possible, to
3852: improve numerical efficiency of matrix-vector products and solves. We
3853: search for consecutive rows with the same nonzero structure, thereby
3854: reusing matrix information to achieve increased efficiency.
3856: Options Database Keys:
3857: + -mat_no_inode - Do not use inodes
3858: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3860: Level: intermediate
3862: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo(),
3863: MatSeqAIJSetTotalPreallocation()
3865: @*/
3866: PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[])
3867: {
3870: PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));
3871: return 0;
3872: }
3874: PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz)
3875: {
3876: Mat_SeqAIJ *b;
3877: PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;
3878: PetscInt i;
3880: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3881: if (nz == MAT_SKIP_ALLOCATION) {
3882: skipallocation = PETSC_TRUE;
3883: nz = 0;
3884: }
3885: PetscLayoutSetUp(B->rmap);
3886: PetscLayoutSetUp(B->cmap);
3888: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3890: if (PetscUnlikelyDebug(nnz)) {
3891: for (i=0; i<B->rmap->n; i++) {
3894: }
3895: }
3897: B->preallocated = PETSC_TRUE;
3899: b = (Mat_SeqAIJ*)B->data;
3901: if (!skipallocation) {
3902: if (!b->imax) {
3903: PetscMalloc1(B->rmap->n,&b->imax);
3904: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3905: }
3906: if (!b->ilen) {
3907: /* b->ilen will count nonzeros in each row so far. */
3908: PetscCalloc1(B->rmap->n,&b->ilen);
3909: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3910: } else {
3911: PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt));
3912: }
3913: if (!b->ipre) {
3914: PetscMalloc1(B->rmap->n,&b->ipre);
3915: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
3916: }
3917: if (!nnz) {
3918: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10;
3919: else if (nz < 0) nz = 1;
3920: nz = PetscMin(nz,B->cmap->n);
3921: for (i=0; i<B->rmap->n; i++) b->imax[i] = nz;
3922: nz = nz*B->rmap->n;
3923: } else {
3924: PetscInt64 nz64 = 0;
3925: for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
3926: PetscIntCast(nz64,&nz);
3927: }
3929: /* allocate the matrix space */
3930: /* FIXME: should B's old memory be unlogged? */
3931: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
3932: if (B->structure_only) {
3933: PetscMalloc1(nz,&b->j);
3934: PetscMalloc1(B->rmap->n+1,&b->i);
3935: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));
3936: } else {
3937: PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);
3938: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
3939: }
3940: b->i[0] = 0;
3941: for (i=1; i<B->rmap->n+1; i++) {
3942: b->i[i] = b->i[i-1] + b->imax[i-1];
3943: }
3944: if (B->structure_only) {
3945: b->singlemalloc = PETSC_FALSE;
3946: b->free_a = PETSC_FALSE;
3947: } else {
3948: b->singlemalloc = PETSC_TRUE;
3949: b->free_a = PETSC_TRUE;
3950: }
3951: b->free_ij = PETSC_TRUE;
3952: } else {
3953: b->free_a = PETSC_FALSE;
3954: b->free_ij = PETSC_FALSE;
3955: }
3957: if (b->ipre && nnz != b->ipre && b->imax) {
3958: /* reserve user-requested sparsity */
3959: PetscArraycpy(b->ipre,b->imax,B->rmap->n);
3960: }
3962: b->nz = 0;
3963: b->maxnz = nz;
3964: B->info.nz_unneeded = (double)b->maxnz;
3965: if (realalloc) {
3966: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
3967: }
3968: B->was_assembled = PETSC_FALSE;
3969: B->assembled = PETSC_FALSE;
3970: return 0;
3971: }
3973: PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A)
3974: {
3975: Mat_SeqAIJ *a;
3976: PetscInt i;
3980: /* Check local size. If zero, then return */
3981: if (!A->rmap->n) return 0;
3983: a = (Mat_SeqAIJ*)A->data;
3984: /* if no saved info, we error out */
3989: PetscArraycpy(a->imax,a->ipre,A->rmap->n);
3990: PetscArrayzero(a->ilen,A->rmap->n);
3991: a->i[0] = 0;
3992: for (i=1; i<A->rmap->n+1; i++) {
3993: a->i[i] = a->i[i-1] + a->imax[i-1];
3994: }
3995: A->preallocated = PETSC_TRUE;
3996: a->nz = 0;
3997: a->maxnz = a->i[A->rmap->n];
3998: A->info.nz_unneeded = (double)a->maxnz;
3999: A->was_assembled = PETSC_FALSE;
4000: A->assembled = PETSC_FALSE;
4001: return 0;
4002: }
4004: /*@
4005: MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format.
4007: Input Parameters:
4008: + B - the matrix
4009: . i - the indices into j for the start of each row (starts with zero)
4010: . j - the column indices for each row (starts with zero) these must be sorted for each row
4011: - v - optional values in the matrix
4013: Level: developer
4015: Notes:
4016: The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays()
4018: This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero
4019: structure will be the union of all the previous nonzero structures.
4021: Developer Notes:
4022: 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
4023: then just copies the v values directly with PetscMemcpy().
4025: This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them.
4027: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation()
4028: @*/
4029: PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[])
4030: {
4033: PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));
4034: return 0;
4035: }
4037: PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[])
4038: {
4039: PetscInt i;
4040: PetscInt m,n;
4041: PetscInt nz;
4042: PetscInt *nnz;
4046: PetscLayoutSetUp(B->rmap);
4047: PetscLayoutSetUp(B->cmap);
4049: MatGetSize(B, &m, &n);
4050: PetscMalloc1(m+1, &nnz);
4051: for (i = 0; i < m; i++) {
4052: nz = Ii[i+1]- Ii[i];
4054: nnz[i] = nz;
4055: }
4056: MatSeqAIJSetPreallocation(B, 0, nnz);
4057: PetscFree(nnz);
4059: for (i = 0; i < m; i++) {
4060: MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES);
4061: }
4063: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
4064: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
4066: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
4067: return 0;
4068: }
4070: /*@
4071: MatSeqAIJKron - Computes C, the Kronecker product of A and B.
4073: Input Parameters:
4074: + A - left-hand side matrix
4075: . B - right-hand side matrix
4076: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4078: Output Parameter:
4079: . C - Kronecker product of A and B
4081: Level: intermediate
4083: Notes:
4084: MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron().
4086: .seealso: MatCreateSeqAIJ(), MATSEQAIJ, MATKAIJ, MatReuse
4087: @*/
4088: PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C)
4089: {
4095: if (reuse == MAT_REUSE_MATRIX) {
4098: }
4099: PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C));
4100: return 0;
4101: }
4103: PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C)
4104: {
4105: Mat newmat;
4106: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4107: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
4108: PetscScalar *v;
4109: const PetscScalar *aa,*ba;
4110: 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;
4111: PetscBool flg;
4117: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg);
4120: if (reuse == MAT_INITIAL_MATRIX) {
4121: PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j);
4122: MatCreate(PETSC_COMM_SELF,&newmat);
4123: MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn);
4124: MatSetType(newmat,MATAIJ);
4125: i[0] = 0;
4126: for (m = 0; m < am; ++m) {
4127: for (p = 0; p < bm; ++p) {
4128: i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]);
4129: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4130: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4131: j[nnz++] = a->j[n]*bn + b->j[q];
4132: }
4133: }
4134: }
4135: }
4136: MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL);
4137: *C = newmat;
4138: PetscFree2(i,j);
4139: nnz = 0;
4140: }
4141: MatSeqAIJGetArray(*C,&v);
4142: MatSeqAIJGetArrayRead(A,&aa);
4143: MatSeqAIJGetArrayRead(B,&ba);
4144: for (m = 0; m < am; ++m) {
4145: for (p = 0; p < bm; ++p) {
4146: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4147: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4148: v[nnz++] = aa[n] * ba[q];
4149: }
4150: }
4151: }
4152: }
4153: MatSeqAIJRestoreArray(*C,&v);
4154: MatSeqAIJRestoreArrayRead(A,&aa);
4155: MatSeqAIJRestoreArrayRead(B,&ba);
4156: return 0;
4157: }
4159: #include <../src/mat/impls/dense/seq/dense.h>
4160: #include <petsc/private/kernels/petscaxpy.h>
4162: /*
4163: Computes (B'*A')' since computing B*A directly is untenable
4165: n p p
4166: [ ] [ ] [ ]
4167: m [ A ] * n [ B ] = m [ C ]
4168: [ ] [ ] [ ]
4170: */
4171: PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C)
4172: {
4173: Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data;
4174: Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data;
4175: Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data;
4176: PetscInt i,j,n,m,q,p;
4177: const PetscInt *ii,*idx;
4178: const PetscScalar *b,*a,*a_q;
4179: PetscScalar *c,*c_q;
4180: PetscInt clda = sub_c->lda;
4181: PetscInt alda = sub_a->lda;
4183: m = A->rmap->n;
4184: n = A->cmap->n;
4185: p = B->cmap->n;
4186: a = sub_a->v;
4187: b = sub_b->a;
4188: c = sub_c->v;
4189: if (clda == m) {
4190: PetscArrayzero(c,m*p);
4191: } else {
4192: for (j=0;j<p;j++)
4193: for (i=0;i<m;i++)
4194: c[j*clda + i] = 0.0;
4195: }
4196: ii = sub_b->i;
4197: idx = sub_b->j;
4198: for (i=0; i<n; i++) {
4199: q = ii[i+1] - ii[i];
4200: while (q-->0) {
4201: c_q = c + clda*(*idx);
4202: a_q = a + alda*i;
4203: PetscKernelAXPY(c_q,*b,a_q,m);
4204: idx++;
4205: b++;
4206: }
4207: }
4208: return 0;
4209: }
4211: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
4212: {
4213: PetscInt m=A->rmap->n,n=B->cmap->n;
4214: PetscBool cisdense;
4217: MatSetSizes(C,m,n,m,n);
4218: MatSetBlockSizesFromMats(C,A,B);
4219: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
4220: if (!cisdense) {
4221: MatSetType(C,MATDENSE);
4222: }
4223: MatSetUp(C);
4225: C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ;
4226: return 0;
4227: }
4229: /* ----------------------------------------------------------------*/
4230: /*MC
4231: MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices,
4232: based on compressed sparse row format.
4234: Options Database Keys:
4235: . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions()
4237: Level: beginner
4239: Notes:
4240: MatSetValues() may be called for this matrix type with a NULL argument for the numerical values,
4241: in this case the values associated with the rows and columns one passes in are set to zero
4242: in the matrix
4244: MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no
4245: space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored
4247: Developer Notes:
4248: It would be nice if all matrix formats supported passing NULL in for the numerical values
4250: .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType, MATSELL, MATSEQSELL, MATMPISELL
4251: M*/
4253: /*MC
4254: MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices.
4256: This matrix type is identical to MATSEQAIJ when constructed with a single process communicator,
4257: and MATMPIAIJ otherwise. As a result, for single process communicators,
4258: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
4259: for communicators controlling multiple processes. It is recommended that you call both of
4260: the above preallocation routines for simplicity.
4262: Options Database Keys:
4263: . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions()
4265: Developer Notes:
4266: Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when
4267: enough exist.
4269: Level: beginner
4271: .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ, MATMPIAIJ, MATSELL, MATSEQSELL, MATMPISELL
4272: M*/
4274: /*MC
4275: MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices.
4277: This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator,
4278: and MATMPIAIJCRL otherwise. As a result, for single process communicators,
4279: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
4280: for communicators controlling multiple processes. It is recommended that you call both of
4281: the above preallocation routines for simplicity.
4283: Options Database Keys:
4284: . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions()
4286: Level: beginner
4288: .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL
4289: M*/
4291: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*);
4292: #if defined(PETSC_HAVE_ELEMENTAL)
4293: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
4294: #endif
4295: #if defined(PETSC_HAVE_SCALAPACK)
4296: PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
4297: #endif
4298: #if defined(PETSC_HAVE_HYPRE)
4299: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*);
4300: #endif
4302: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*);
4303: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*);
4304: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat);
4306: /*@C
4307: MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored
4309: Not Collective
4311: Input Parameter:
4312: . mat - a MATSEQAIJ matrix
4314: Output Parameter:
4315: . array - pointer to the data
4317: Level: intermediate
4319: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4320: @*/
4321: PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array)
4322: {
4323: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4325: if (aij->ops->getarray) {
4326: (*aij->ops->getarray)(A,array);
4327: } else {
4328: *array = aij->a;
4329: }
4330: return 0;
4331: }
4333: /*@C
4334: MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray()
4336: Not Collective
4338: Input Parameters:
4339: + mat - a MATSEQAIJ matrix
4340: - array - pointer to the data
4342: Level: intermediate
4344: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90()
4345: @*/
4346: PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array)
4347: {
4348: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4350: if (aij->ops->restorearray) {
4351: (*aij->ops->restorearray)(A,array);
4352: } else {
4353: *array = NULL;
4354: }
4355: MatSeqAIJInvalidateDiagonal(A);
4356: PetscObjectStateIncrease((PetscObject)A);
4357: return 0;
4358: }
4360: /*@C
4361: MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored
4363: Not Collective
4365: Input Parameter:
4366: . mat - a MATSEQAIJ matrix
4368: Output Parameter:
4369: . array - pointer to the data
4371: Level: intermediate
4373: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead()
4374: @*/
4375: PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array)
4376: {
4377: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4379: if (aij->ops->getarrayread) {
4380: (*aij->ops->getarrayread)(A,array);
4381: } else {
4382: *array = aij->a;
4383: }
4384: return 0;
4385: }
4387: /*@C
4388: MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead
4390: Not Collective
4392: Input Parameter:
4393: . mat - a MATSEQAIJ matrix
4395: Output Parameter:
4396: . array - pointer to the data
4398: Level: intermediate
4400: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4401: @*/
4402: PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array)
4403: {
4404: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4406: if (aij->ops->restorearrayread) {
4407: (*aij->ops->restorearrayread)(A,array);
4408: } else {
4409: *array = NULL;
4410: }
4411: return 0;
4412: }
4414: /*@C
4415: MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored
4417: Not Collective
4419: Input Parameter:
4420: . mat - a MATSEQAIJ matrix
4422: Output Parameter:
4423: . array - pointer to the data
4425: Level: intermediate
4427: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead()
4428: @*/
4429: PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array)
4430: {
4431: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4433: if (aij->ops->getarraywrite) {
4434: (*aij->ops->getarraywrite)(A,array);
4435: } else {
4436: *array = aij->a;
4437: }
4438: MatSeqAIJInvalidateDiagonal(A);
4439: PetscObjectStateIncrease((PetscObject)A);
4440: return 0;
4441: }
4443: /*@C
4444: MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead
4446: Not Collective
4448: Input Parameter:
4449: . mat - a MATSEQAIJ matrix
4451: Output Parameter:
4452: . array - pointer to the data
4454: Level: intermediate
4456: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4457: @*/
4458: PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array)
4459: {
4460: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4462: if (aij->ops->restorearraywrite) {
4463: (*aij->ops->restorearraywrite)(A,array);
4464: } else {
4465: *array = NULL;
4466: }
4467: return 0;
4468: }
4470: /*@C
4471: MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix
4473: Not Collective
4475: Input Parameter:
4476: . mat - a matrix of type MATSEQAIJ or its subclasses
4478: Output Parameters:
4479: + i - row map array of the matrix
4480: . j - column index array of the matrix
4481: . a - data array of the matrix
4482: - memtype - memory type of the arrays
4484: Notes:
4485: Any of the output parameters can be NULL, in which case the corresponding value is not returned.
4486: If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host.
4488: One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix.
4489: If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix.
4491: Level: Developer
4493: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4494: @*/
4495: PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype)
4496: {
4497: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
4500: if (aij->ops->getcsrandmemtype) {
4501: (*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype);
4502: } else {
4503: if (i) *i = aij->i;
4504: if (j) *j = aij->j;
4505: if (a) *a = aij->a;
4506: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
4507: }
4508: return 0;
4509: }
4511: /*@C
4512: MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row
4514: Not Collective
4516: Input Parameter:
4517: . mat - a MATSEQAIJ matrix
4519: Output Parameter:
4520: . nz - the maximum number of nonzeros in any row
4522: Level: intermediate
4524: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4525: @*/
4526: PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz)
4527: {
4528: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4530: *nz = aij->rmax;
4531: return 0;
4532: }
4534: PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[])
4535: {
4536: MPI_Comm comm;
4537: PetscInt *i,*j;
4538: PetscInt M,N,row;
4539: PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */
4540: PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */
4541: PetscInt *Aj;
4542: PetscScalar *Aa;
4543: Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data);
4544: MatType rtype;
4545: PetscCount *perm,*jmap;
4547: MatResetPreallocationCOO_SeqAIJ(mat);
4548: PetscObjectGetComm((PetscObject)mat,&comm);
4549: MatGetSize(mat,&M,&N);
4550: PetscMalloc2(coo_n,&i,coo_n,&j);
4551: PetscArraycpy(i,coo_i,coo_n); /* Make a copy since we'll modify it */
4552: PetscArraycpy(j,coo_j,coo_n);
4553: PetscMalloc1(coo_n,&perm);
4554: for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */
4555: if (j[k] < 0) i[k] = -1;
4556: perm[k] = k;
4557: }
4559: /* Sort by row */
4560: PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm);
4561: for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */
4562: nneg = k;
4563: 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 */
4564: nnz = 0; /* Total number of unique nonzeros to be counted */
4565: jmap++; /* Inc jmap by 1 for convinience */
4567: PetscCalloc1(M+1,&Ai); /* CSR of A */
4568: PetscMalloc1(coo_n-nneg,&Aj); /* We have at most coo_n-nneg unique nonzeros */
4570: /* In each row, sort by column, then unique column indices to get row length */
4571: Ai++; /* Inc by 1 for convinience */
4572: q = 0; /* q-th unique nonzero, with q starting from 0 */
4573: while (k<coo_n) {
4574: row = i[k];
4575: start = k; /* [start,end) indices for this row */
4576: while (k<coo_n && i[k] == row) k++;
4577: end = k;
4578: PetscSortIntWithCountArray(end-start,j+start,perm+start);
4579: /* Find number of unique col entries in this row */
4580: Aj[q] = j[start]; /* Log the first nonzero in this row */
4581: jmap[q] = 1; /* Number of repeats of this nozero entry */
4582: Ai[row] = 1;
4583: nnz++;
4585: for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */
4586: if (j[p] != j[p-1]) { /* Meet a new nonzero */
4587: q++;
4588: jmap[q] = 1;
4589: Aj[q] = j[p];
4590: Ai[row]++;
4591: nnz++;
4592: } else {
4593: jmap[q]++;
4594: }
4595: }
4596: q++; /* Move to next row and thus next unique nonzero */
4597: }
4598: PetscFree2(i,j);
4600: Ai--; /* Back to the beginning of Ai[] */
4601: for (k=0; k<M; k++) Ai[k+1] += Ai[k];
4602: jmap--; /* Back to the beginning of jmap[] */
4603: jmap[0] = 0;
4604: for (k=0; k<nnz; k++) jmap[k+1] += jmap[k];
4605: if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */
4606: PetscCount *jmap_new;
4607: PetscInt *Aj_new;
4609: PetscMalloc1(nnz+1,&jmap_new);
4610: PetscArraycpy(jmap_new,jmap,nnz+1);
4611: PetscFree(jmap);
4612: jmap = jmap_new;
4614: PetscMalloc1(nnz,&Aj_new);
4615: PetscArraycpy(Aj_new,Aj,nnz);
4616: PetscFree(Aj);
4617: Aj = Aj_new;
4618: }
4620: if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */
4621: PetscCount *perm_new;
4623: PetscMalloc1(coo_n-nneg,&perm_new);
4624: PetscArraycpy(perm_new,perm+nneg,coo_n-nneg);
4625: PetscFree(perm);
4626: perm = perm_new;
4627: }
4629: MatGetRootType_Private(mat,&rtype);
4630: PetscCalloc1(nnz,&Aa); /* Zero the matrix */
4631: MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat);
4633: seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */
4634: seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */
4635: /* Record COO fields */
4636: seqaij->coo_n = coo_n;
4637: seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */
4638: seqaij->jmap = jmap; /* of length nnz+1 */
4639: seqaij->perm = perm;
4640: return 0;
4641: }
4643: static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode)
4644: {
4645: Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data;
4646: PetscCount i,j,Annz = aseq->nz;
4647: PetscCount *perm = aseq->perm,*jmap = aseq->jmap;
4648: PetscScalar *Aa;
4650: MatSeqAIJGetArray(A,&Aa);
4651: for (i=0; i<Annz; i++) {
4652: PetscScalar sum = 0.0;
4653: for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]];
4654: Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum;
4655: }
4656: MatSeqAIJRestoreArray(A,&Aa);
4657: return 0;
4658: }
4660: #if defined(PETSC_HAVE_CUDA)
4661: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*);
4662: #endif
4663: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4664: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*);
4665: #endif
4667: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B)
4668: {
4669: Mat_SeqAIJ *b;
4670: PetscMPIInt size;
4672: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
4675: PetscNewLog(B,&b);
4677: B->data = (void*)b;
4679: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
4680: if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull;
4682: b->row = NULL;
4683: b->col = NULL;
4684: b->icol = NULL;
4685: b->reallocs = 0;
4686: b->ignorezeroentries = PETSC_FALSE;
4687: b->roworiented = PETSC_TRUE;
4688: b->nonew = 0;
4689: b->diag = NULL;
4690: b->solve_work = NULL;
4691: B->spptr = NULL;
4692: b->saved_values = NULL;
4693: b->idiag = NULL;
4694: b->mdiag = NULL;
4695: b->ssor_work = NULL;
4696: b->omega = 1.0;
4697: b->fshift = 0.0;
4698: b->idiagvalid = PETSC_FALSE;
4699: b->ibdiagvalid = PETSC_FALSE;
4700: b->keepnonzeropattern = PETSC_FALSE;
4702: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4704: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4705: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);
4706: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);
4707: #endif
4709: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);
4710: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);
4711: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);
4712: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);
4713: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);
4714: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);
4715: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);
4716: #if defined(PETSC_HAVE_MKL_SPARSE)
4717: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);
4718: #endif
4719: #if defined(PETSC_HAVE_CUDA)
4720: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);
4721: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4722: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ);
4723: #endif
4724: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4725: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos);
4726: #endif
4727: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);
4728: #if defined(PETSC_HAVE_ELEMENTAL)
4729: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);
4730: #endif
4731: #if defined(PETSC_HAVE_SCALAPACK)
4732: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK);
4733: #endif
4734: #if defined(PETSC_HAVE_HYPRE)
4735: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);
4736: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ);
4737: #endif
4738: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);
4739: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);
4740: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);
4741: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);
4742: PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);
4743: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);
4744: PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);
4745: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);
4746: PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);
4747: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ);
4748: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ);
4749: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4750: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ);
4751: PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ);
4752: PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ);
4753: MatCreate_SeqAIJ_Inode(B);
4754: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4755: MatSeqAIJSetTypeFromOptions(B); /* this allows changing the matrix subtype to say MATSEQAIJPERM */
4756: return 0;
4757: }
4759: /*
4760: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
4761: */
4762: PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
4763: {
4764: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data;
4765: PetscInt m = A->rmap->n,i;
4769: C->factortype = A->factortype;
4770: c->row = NULL;
4771: c->col = NULL;
4772: c->icol = NULL;
4773: c->reallocs = 0;
4775: C->assembled = A->assembled;
4776: C->preallocated = A->preallocated;
4778: if (A->preallocated) {
4779: PetscLayoutReference(A->rmap,&C->rmap);
4780: PetscLayoutReference(A->cmap,&C->cmap);
4782: PetscMalloc1(m,&c->imax);
4783: PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt));
4784: PetscMalloc1(m,&c->ilen);
4785: PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt));
4786: PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));
4788: /* allocate the matrix space */
4789: if (mallocmatspace) {
4790: PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);
4791: PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));
4793: c->singlemalloc = PETSC_TRUE;
4795: PetscArraycpy(c->i,a->i,m+1);
4796: if (m > 0) {
4797: PetscArraycpy(c->j,a->j,a->i[m]);
4798: if (cpvalues == MAT_COPY_VALUES) {
4799: const PetscScalar *aa;
4801: MatSeqAIJGetArrayRead(A,&aa);
4802: PetscArraycpy(c->a,aa,a->i[m]);
4803: MatSeqAIJGetArrayRead(A,&aa);
4804: } else {
4805: PetscArrayzero(c->a,a->i[m]);
4806: }
4807: }
4808: }
4810: c->ignorezeroentries = a->ignorezeroentries;
4811: c->roworiented = a->roworiented;
4812: c->nonew = a->nonew;
4813: if (a->diag) {
4814: PetscMalloc1(m+1,&c->diag);
4815: PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt));
4816: PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));
4817: } else c->diag = NULL;
4819: c->solve_work = NULL;
4820: c->saved_values = NULL;
4821: c->idiag = NULL;
4822: c->ssor_work = NULL;
4823: c->keepnonzeropattern = a->keepnonzeropattern;
4824: c->free_a = PETSC_TRUE;
4825: c->free_ij = PETSC_TRUE;
4827: c->rmax = a->rmax;
4828: c->nz = a->nz;
4829: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
4831: c->compressedrow.use = a->compressedrow.use;
4832: c->compressedrow.nrows = a->compressedrow.nrows;
4833: if (a->compressedrow.use) {
4834: i = a->compressedrow.nrows;
4835: PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);
4836: PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);
4837: PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);
4838: } else {
4839: c->compressedrow.use = PETSC_FALSE;
4840: c->compressedrow.i = NULL;
4841: c->compressedrow.rindex = NULL;
4842: }
4843: c->nonzerorowcnt = a->nonzerorowcnt;
4844: C->nonzerostate = A->nonzerostate;
4846: MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);
4847: }
4848: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
4849: return 0;
4850: }
4852: PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
4853: {
4854: MatCreate(PetscObjectComm((PetscObject)A),B);
4855: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
4856: if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) {
4857: MatSetBlockSizesFromMats(*B,A,A);
4858: }
4859: MatSetType(*B,((PetscObject)A)->type_name);
4860: MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);
4861: return 0;
4862: }
4864: PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer)
4865: {
4866: PetscBool isbinary, ishdf5;
4870: /* force binary viewer to load .info file if it has not yet done so */
4871: PetscViewerSetUp(viewer);
4872: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
4873: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
4874: if (isbinary) {
4875: MatLoad_SeqAIJ_Binary(newMat,viewer);
4876: } else if (ishdf5) {
4877: #if defined(PETSC_HAVE_HDF5)
4878: MatLoad_AIJ_HDF5(newMat,viewer);
4879: #else
4880: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
4881: #endif
4882: } else {
4883: 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);
4884: }
4885: return 0;
4886: }
4888: PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer)
4889: {
4890: Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data;
4891: PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i;
4893: PetscViewerSetUp(viewer);
4895: /* read in matrix header */
4896: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
4898: M = header[1]; N = header[2]; nz = header[3];
4903: /* set block sizes from the viewer's .info file */
4904: MatLoad_Binary_BlockSizes(mat,viewer);
4905: /* set local and global sizes if not set already */
4906: if (mat->rmap->n < 0) mat->rmap->n = M;
4907: if (mat->cmap->n < 0) mat->cmap->n = N;
4908: if (mat->rmap->N < 0) mat->rmap->N = M;
4909: if (mat->cmap->N < 0) mat->cmap->N = N;
4910: PetscLayoutSetUp(mat->rmap);
4911: PetscLayoutSetUp(mat->cmap);
4913: /* check if the matrix sizes are correct */
4914: MatGetSize(mat,&rows,&cols);
4917: /* read in row lengths */
4918: PetscMalloc1(M,&rowlens);
4919: PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT);
4920: /* check if sum(rowlens) is same as nz */
4921: sum = 0; for (i=0; i<M; i++) sum += rowlens[i];
4923: /* preallocate and check sizes */
4924: MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens);
4925: MatGetSize(mat,&rows,&cols);
4927: /* store row lengths */
4928: PetscArraycpy(a->ilen,rowlens,M);
4929: PetscFree(rowlens);
4931: /* fill in "i" row pointers */
4932: a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i];
4933: /* read in "j" column indices */
4934: PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT);
4935: /* read in "a" nonzero values */
4936: PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR);
4938: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
4939: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
4940: return 0;
4941: }
4943: PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg)
4944: {
4945: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data;
4946: const PetscScalar *aa,*ba;
4947: #if defined(PETSC_USE_COMPLEX)
4948: PetscInt k;
4949: #endif
4951: /* If the matrix dimensions are not equal,or no of nonzeros */
4952: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) {
4953: *flg = PETSC_FALSE;
4954: return 0;
4955: }
4957: /* if the a->i are the same */
4958: PetscArraycmp(a->i,b->i,A->rmap->n+1,flg);
4959: if (!*flg) return 0;
4961: /* if a->j are the same */
4962: PetscArraycmp(a->j,b->j,a->nz,flg);
4963: if (!*flg) return 0;
4965: MatSeqAIJGetArrayRead(A,&aa);
4966: MatSeqAIJGetArrayRead(B,&ba);
4967: /* if a->a are the same */
4968: #if defined(PETSC_USE_COMPLEX)
4969: for (k=0; k<a->nz; k++) {
4970: if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) {
4971: *flg = PETSC_FALSE;
4972: return 0;
4973: }
4974: }
4975: #else
4976: PetscArraycmp(aa,ba,a->nz,flg);
4977: #endif
4978: MatSeqAIJRestoreArrayRead(A,&aa);
4979: MatSeqAIJRestoreArrayRead(B,&ba);
4980: return 0;
4981: }
4983: /*@
4984: MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format)
4985: provided by the user.
4987: Collective
4989: Input Parameters:
4990: + comm - must be an MPI communicator of size 1
4991: . m - number of rows
4992: . n - number of columns
4993: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix
4994: . j - column indices
4995: - a - matrix values
4997: Output Parameter:
4998: . mat - the matrix
5000: Level: intermediate
5002: Notes:
5003: The i, j, and a arrays are not copied by this routine, the user must free these arrays
5004: once the matrix is destroyed and not before
5006: You cannot set new nonzero locations into this matrix, that will generate an error.
5008: The i and j indices are 0 based
5010: The format which is used for the sparse matrix input, is equivalent to a
5011: row-major ordering.. i.e for the following matrix, the input data expected is
5012: as shown
5014: $ 1 0 0
5015: $ 2 0 3
5016: $ 4 5 6
5017: $
5018: $ i = {0,1,3,6} [size = nrow+1 = 3+1]
5019: $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row
5020: $ v = {1,2,3,4,5,6} [size = 6]
5022: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5024: @*/
5025: PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
5026: {
5027: PetscInt ii;
5028: Mat_SeqAIJ *aij;
5029: PetscInt jj;
5032: MatCreate(comm,mat);
5033: MatSetSizes(*mat,m,n,m,n);
5034: /* MatSetBlockSizes(*mat,,); */
5035: MatSetType(*mat,MATSEQAIJ);
5036: MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL);
5037: aij = (Mat_SeqAIJ*)(*mat)->data;
5038: PetscMalloc1(m,&aij->imax);
5039: PetscMalloc1(m,&aij->ilen);
5041: aij->i = i;
5042: aij->j = j;
5043: aij->a = a;
5044: aij->singlemalloc = PETSC_FALSE;
5045: aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
5046: aij->free_a = PETSC_FALSE;
5047: aij->free_ij = PETSC_FALSE;
5049: for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) {
5050: aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
5051: if (PetscDefined(USE_DEBUG)) {
5053: for (jj=i[ii]+1; jj<i[ii+1]; jj++) {
5056: }
5057: }
5058: }
5059: if (PetscDefined(USE_DEBUG)) {
5060: for (ii=0; ii<aij->i[m]; ii++) {
5063: }
5064: }
5066: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5067: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5068: return 0;
5069: }
5071: /*@C
5072: MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format)
5073: provided by the user.
5075: Collective
5077: Input Parameters:
5078: + comm - must be an MPI communicator of size 1
5079: . m - number of rows
5080: . n - number of columns
5081: . i - row indices
5082: . j - column indices
5083: . a - matrix values
5084: . nz - number of nonzeros
5085: - idx - 0 or 1 based
5087: Output Parameter:
5088: . mat - the matrix
5090: Level: intermediate
5092: Notes:
5093: 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,
5094: the input data expected is as shown
5095: .vb
5096: 1 0 0
5097: 2 0 3
5098: 4 5 6
5100: i = {0,1,1,2,2,2}
5101: j = {0,0,2,0,1,2}
5102: v = {1,2,3,4,5,6}
5103: .ve
5105: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5107: @*/
5108: PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx)
5109: {
5110: PetscInt ii, *nnz, one = 1,row,col;
5112: PetscCalloc1(m,&nnz);
5113: for (ii = 0; ii < nz; ii++) {
5114: nnz[i[ii] - !!idx] += 1;
5115: }
5116: MatCreate(comm,mat);
5117: MatSetSizes(*mat,m,n,m,n);
5118: MatSetType(*mat,MATSEQAIJ);
5119: MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);
5120: for (ii = 0; ii < nz; ii++) {
5121: if (idx) {
5122: row = i[ii] - 1;
5123: col = j[ii] - 1;
5124: } else {
5125: row = i[ii];
5126: col = j[ii];
5127: }
5128: MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);
5129: }
5130: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5131: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5132: PetscFree(nnz);
5133: return 0;
5134: }
5136: PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A)
5137: {
5138: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
5140: a->idiagvalid = PETSC_FALSE;
5141: a->ibdiagvalid = PETSC_FALSE;
5143: MatSeqAIJInvalidateDiagonal_Inode(A);
5144: return 0;
5145: }
5147: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
5148: {
5149: PetscMPIInt size;
5151: MPI_Comm_size(comm,&size);
5152: if (size == 1) {
5153: if (scall == MAT_INITIAL_MATRIX) {
5154: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
5155: } else {
5156: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
5157: }
5158: } else {
5159: MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);
5160: }
5161: return 0;
5162: }
5164: /*
5165: Permute A into C's *local* index space using rowemb,colemb.
5166: The embedding are supposed to be injections and the above implies that the range of rowemb is a subset
5167: of [0,m), colemb is in [0,n).
5168: If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A.
5169: */
5170: PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B)
5171: {
5172: /* If making this function public, change the error returned in this function away from _PLIB. */
5173: Mat_SeqAIJ *Baij;
5174: PetscBool seqaij;
5175: PetscInt m,n,*nz,i,j,count;
5176: PetscScalar v;
5177: const PetscInt *rowindices,*colindices;
5179: if (!B) return 0;
5180: /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */
5181: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
5183: if (rowemb) {
5184: ISGetLocalSize(rowemb,&m);
5186: } else {
5188: }
5189: if (colemb) {
5190: ISGetLocalSize(colemb,&n);
5192: } else {
5194: }
5196: Baij = (Mat_SeqAIJ*)(B->data);
5197: if (pattern == DIFFERENT_NONZERO_PATTERN) {
5198: PetscMalloc1(B->rmap->n,&nz);
5199: for (i=0; i<B->rmap->n; i++) {
5200: nz[i] = Baij->i[i+1] - Baij->i[i];
5201: }
5202: MatSeqAIJSetPreallocation(C,0,nz);
5203: PetscFree(nz);
5204: }
5205: if (pattern == SUBSET_NONZERO_PATTERN) {
5206: MatZeroEntries(C);
5207: }
5208: count = 0;
5209: rowindices = NULL;
5210: colindices = NULL;
5211: if (rowemb) {
5212: ISGetIndices(rowemb,&rowindices);
5213: }
5214: if (colemb) {
5215: ISGetIndices(colemb,&colindices);
5216: }
5217: for (i=0; i<B->rmap->n; i++) {
5218: PetscInt row;
5219: row = i;
5220: if (rowindices) row = rowindices[i];
5221: for (j=Baij->i[i]; j<Baij->i[i+1]; j++) {
5222: PetscInt col;
5223: col = Baij->j[count];
5224: if (colindices) col = colindices[col];
5225: v = Baij->a[count];
5226: MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);
5227: ++count;
5228: }
5229: }
5230: /* FIXME: set C's nonzerostate correctly. */
5231: /* Assembly for C is necessary. */
5232: C->preallocated = PETSC_TRUE;
5233: C->assembled = PETSC_TRUE;
5234: C->was_assembled = PETSC_FALSE;
5235: return 0;
5236: }
5238: PetscFunctionList MatSeqAIJList = NULL;
5240: /*@C
5241: MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype
5243: Collective on Mat
5245: Input Parameters:
5246: + mat - the matrix object
5247: - matype - matrix type
5249: Options Database Key:
5250: . -mat_seqai_type <method> - for example seqaijcrl
5252: Level: intermediate
5254: .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat
5255: @*/
5256: PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype)
5257: {
5258: PetscBool sametype;
5259: PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*);
5262: PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);
5263: if (sametype) return 0;
5265: PetscFunctionListFind(MatSeqAIJList,matype,&r);
5267: (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);
5268: return 0;
5269: }
5271: /*@C
5272: MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices
5274: Not Collective
5276: Input Parameters:
5277: + name - name of a new user-defined matrix type, for example MATSEQAIJCRL
5278: - function - routine to convert to subtype
5280: Notes:
5281: MatSeqAIJRegister() may be called multiple times to add several user-defined solvers.
5283: Then, your matrix can be chosen with the procedural interface at runtime via the option
5284: $ -mat_seqaij_type my_mat
5286: Level: advanced
5288: .seealso: MatSeqAIJRegisterAll()
5290: Level: advanced
5291: @*/
5292: PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *))
5293: {
5294: MatInitializePackage();
5295: PetscFunctionListAdd(&MatSeqAIJList,sname,function);
5296: return 0;
5297: }
5299: PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE;
5301: /*@C
5302: MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ
5304: Not Collective
5306: Level: advanced
5308: .seealso: MatRegisterAll(), MatSeqAIJRegister()
5309: @*/
5310: PetscErrorCode MatSeqAIJRegisterAll(void)
5311: {
5312: if (MatSeqAIJRegisterAllCalled) return 0;
5313: MatSeqAIJRegisterAllCalled = PETSC_TRUE;
5315: MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);
5316: MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);
5317: MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);
5318: #if defined(PETSC_HAVE_MKL_SPARSE)
5319: MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);
5320: #endif
5321: #if defined(PETSC_HAVE_CUDA)
5322: MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE);
5323: #endif
5324: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
5325: MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos);
5326: #endif
5327: #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA)
5328: MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);
5329: #endif
5330: return 0;
5331: }
5333: /*
5334: Special version for direct calls from Fortran
5335: */
5336: #include <petsc/private/fortranimpl.h>
5337: #if defined(PETSC_HAVE_FORTRAN_CAPS)
5338: #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ
5339: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
5340: #define matsetvaluesseqaij_ matsetvaluesseqaij
5341: #endif
5343: /* Change these macros so can be used in void function */
5345: /* Change these macros so can be used in void function */
5346: /* Identical to PetscCallVoid, except it assigns to *_ierr */
5347: #undef PetscCall
5348: #define PetscCall(...) do { \
5349: PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \
5350: if (PetscUnlikely(ierr_msv_mpiaij)) { \
5351: *_PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \
5352: return; \
5353: } \
5354: } while (0)
5356: #undef SETERRQ
5357: #define SETERRQ(comm,ierr,...) do { \
5358: *_PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \
5359: return; \
5360: } while (0)
5362: PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr)
5363: {
5364: Mat A = *AA;
5365: PetscInt m = *mm, n = *nn;
5366: InsertMode is = *isis;
5367: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
5368: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
5369: PetscInt *imax,*ai,*ailen;
5370: PetscInt *aj,nonew = a->nonew,lastcol = -1;
5371: MatScalar *ap,value,*aa;
5372: PetscBool ignorezeroentries = a->ignorezeroentries;
5373: PetscBool roworiented = a->roworiented;
5375: MatCheckPreallocated(A,1);
5376: imax = a->imax;
5377: ai = a->i;
5378: ailen = a->ilen;
5379: aj = a->j;
5380: aa = a->a;
5382: for (k=0; k<m; k++) { /* loop over added rows */
5383: row = im[k];
5384: if (row < 0) continue;
5386: rp = aj + ai[row]; ap = aa + ai[row];
5387: rmax = imax[row]; nrow = ailen[row];
5388: low = 0;
5389: high = nrow;
5390: for (l=0; l<n; l++) { /* loop over added columns */
5391: if (in[l] < 0) continue;
5393: col = in[l];
5394: if (roworiented) value = v[l + k*n];
5395: else value = v[k + l*m];
5397: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
5399: if (col <= lastcol) low = 0;
5400: else high = nrow;
5401: lastcol = col;
5402: while (high-low > 5) {
5403: t = (low+high)/2;
5404: if (rp[t] > col) high = t;
5405: else low = t;
5406: }
5407: for (i=low; i<high; i++) {
5408: if (rp[i] > col) break;
5409: if (rp[i] == col) {
5410: if (is == ADD_VALUES) ap[i] += value;
5411: else ap[i] = value;
5412: goto noinsert;
5413: }
5414: }
5415: if (value == 0.0 && ignorezeroentries) goto noinsert;
5416: if (nonew == 1) goto noinsert;
5418: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
5419: N = nrow++ - 1; a->nz++; high++;
5420: /* shift up all the later entries in this row */
5421: for (ii=N; ii>=i; ii--) {
5422: rp[ii+1] = rp[ii];
5423: ap[ii+1] = ap[ii];
5424: }
5425: rp[i] = col;
5426: ap[i] = value;
5427: A->nonzerostate++;
5428: noinsert:;
5429: low = i + 1;
5430: }
5431: ailen[row] = nrow;
5432: }
5433: return;
5434: }
5435: /* Undefining these here since they were redefined from their original definition above! No
5436: * other PETSc functions should be defined past this point, as it is impossible to recover the
5437: * original definitions */
5438: #undef PetscCall
5439: #undef SETERRQ