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