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 && x->nz == y->nz) {
3178: PetscBool e;
3179: PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e);
3180: if (e) {
3181: PetscArraycmp(x->j,y->j,y->nz,&e);
3182: if (e) {
3183: str = SAME_NONZERO_PATTERN;
3184: }
3185: }
3186: }
3187: if (str == SAME_NONZERO_PATTERN) {
3188: const PetscScalar *xa;
3189: PetscScalar *ya,alpha = a;
3190: PetscBLASInt one = 1,bnz;
3192: PetscBLASIntCast(x->nz,&bnz);
3193: MatSeqAIJGetArray(Y,&ya);
3194: MatSeqAIJGetArrayRead(X,&xa);
3195: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one));
3196: MatSeqAIJRestoreArrayRead(X,&xa);
3197: MatSeqAIJRestoreArray(Y,&ya);
3198: PetscLogFlops(2.0*bnz);
3199: MatSeqAIJInvalidateDiagonal(Y);
3200: PetscObjectStateIncrease((PetscObject)Y);
3201: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
3202: MatAXPY_Basic(Y,a,X,str);
3203: } else {
3204: Mat B;
3205: PetscInt *nnz;
3206: PetscMalloc1(Y->rmap->N,&nnz);
3207: MatCreate(PetscObjectComm((PetscObject)Y),&B);
3208: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
3209: MatSetLayouts(B,Y->rmap,Y->cmap);
3210: MatSetType(B,((PetscObject)Y)->type_name);
3211: MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);
3212: MatSeqAIJSetPreallocation(B,0,nnz);
3213: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
3214: MatHeaderReplace(Y,&B);
3215: PetscFree(nnz);
3216: }
3217: return(0);
3218: }
3220: PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat)
3221: {
3222: #if defined(PETSC_USE_COMPLEX)
3223: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3224: PetscInt i,nz;
3225: PetscScalar *a;
3229: nz = aij->nz;
3230: MatSeqAIJGetArray(mat,&a);
3231: for (i=0; i<nz; i++) a[i] = PetscConj(a[i]);
3232: MatSeqAIJRestoreArray(mat,&a);
3233: #else
3235: #endif
3236: return(0);
3237: }
3239: PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3240: {
3241: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3242: PetscErrorCode ierr;
3243: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3244: PetscReal atmp;
3245: PetscScalar *x;
3246: const MatScalar *aa,*av;
3249: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3250: MatSeqAIJGetArrayRead(A,&av);
3251: aa = av;
3252: ai = a->i;
3253: aj = a->j;
3255: VecSet(v,0.0);
3256: VecGetArrayWrite(v,&x);
3257: VecGetLocalSize(v,&n);
3258: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
3259: for (i=0; i<m; i++) {
3260: ncols = ai[1] - ai[0]; ai++;
3261: for (j=0; j<ncols; j++) {
3262: atmp = PetscAbsScalar(*aa);
3263: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
3264: aa++; aj++;
3265: }
3266: }
3267: VecRestoreArrayWrite(v,&x);
3268: MatSeqAIJRestoreArrayRead(A,&av);
3269: return(0);
3270: }
3272: PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3273: {
3274: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3275: PetscErrorCode ierr;
3276: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3277: PetscScalar *x;
3278: const MatScalar *aa,*av;
3281: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3282: MatSeqAIJGetArrayRead(A,&av);
3283: aa = av;
3284: ai = a->i;
3285: aj = a->j;
3287: VecSet(v,0.0);
3288: VecGetArrayWrite(v,&x);
3289: VecGetLocalSize(v,&n);
3290: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
3291: for (i=0; i<m; i++) {
3292: ncols = ai[1] - ai[0]; ai++;
3293: if (ncols == A->cmap->n) { /* row is dense */
3294: x[i] = *aa; if (idx) idx[i] = 0;
3295: } else { /* row is sparse so already KNOW maximum is 0.0 or higher */
3296: x[i] = 0.0;
3297: if (idx) {
3298: for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */
3299: if (aj[j] > j) {
3300: idx[i] = j;
3301: break;
3302: }
3303: }
3304: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3305: if (j==ncols && j < A->cmap->n) idx[i] = j;
3306: }
3307: }
3308: for (j=0; j<ncols; j++) {
3309: if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3310: aa++; aj++;
3311: }
3312: }
3313: VecRestoreArrayWrite(v,&x);
3314: MatSeqAIJRestoreArrayRead(A,&av);
3315: return(0);
3316: }
3318: PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3319: {
3320: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3321: PetscErrorCode ierr;
3322: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
3323: PetscScalar *x;
3324: const MatScalar *aa,*av;
3327: MatSeqAIJGetArrayRead(A,&av);
3328: aa = av;
3329: ai = a->i;
3330: aj = a->j;
3332: VecSet(v,0.0);
3333: VecGetArrayWrite(v,&x);
3334: VecGetLocalSize(v,&n);
3335: if (n != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %D vs. %D rows", m, n);
3336: for (i=0; i<m; i++) {
3337: ncols = ai[1] - ai[0]; ai++;
3338: if (ncols == A->cmap->n) { /* row is dense */
3339: x[i] = *aa; if (idx) idx[i] = 0;
3340: } else { /* row is sparse so already KNOW minimum is 0.0 or higher */
3341: x[i] = 0.0;
3342: if (idx) { /* find first implicit 0.0 in the row */
3343: for (j=0; j<ncols; j++) {
3344: if (aj[j] > j) {
3345: idx[i] = j;
3346: break;
3347: }
3348: }
3349: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3350: if (j==ncols && j < A->cmap->n) idx[i] = j;
3351: }
3352: }
3353: for (j=0; j<ncols; j++) {
3354: if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3355: aa++; aj++;
3356: }
3357: }
3358: VecRestoreArrayWrite(v,&x);
3359: MatSeqAIJRestoreArrayRead(A,&av);
3360: return(0);
3361: }
3363: PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[])
3364: {
3365: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3366: PetscErrorCode ierr;
3367: PetscInt i,j,m = A->rmap->n,ncols,n;
3368: const PetscInt *ai,*aj;
3369: PetscScalar *x;
3370: const MatScalar *aa,*av;
3373: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3374: MatSeqAIJGetArrayRead(A,&av);
3375: aa = av;
3376: ai = a->i;
3377: aj = a->j;
3379: VecSet(v,0.0);
3380: VecGetArrayWrite(v,&x);
3381: VecGetLocalSize(v,&n);
3382: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
3383: for (i=0; i<m; i++) {
3384: ncols = ai[1] - ai[0]; ai++;
3385: if (ncols == A->cmap->n) { /* row is dense */
3386: x[i] = *aa; if (idx) idx[i] = 0;
3387: } else { /* row is sparse so already KNOW minimum is 0.0 or lower */
3388: x[i] = 0.0;
3389: if (idx) { /* find first implicit 0.0 in the row */
3390: for (j=0; j<ncols; j++) {
3391: if (aj[j] > j) {
3392: idx[i] = j;
3393: break;
3394: }
3395: }
3396: /* in case first implicit 0.0 in the row occurs at ncols-th column */
3397: if (j==ncols && j < A->cmap->n) idx[i] = j;
3398: }
3399: }
3400: for (j=0; j<ncols; j++) {
3401: if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
3402: aa++; aj++;
3403: }
3404: }
3405: VecRestoreArrayWrite(v,&x);
3406: MatSeqAIJRestoreArrayRead(A,&av);
3407: return(0);
3408: }
3410: PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values)
3411: {
3412: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
3413: PetscErrorCode ierr;
3414: PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j;
3415: MatScalar *diag,work[25],*v_work;
3416: const PetscReal shift = 0.0;
3417: PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE;
3420: allowzeropivot = PetscNot(A->erroriffailure);
3421: if (a->ibdiagvalid) {
3422: if (values) *values = a->ibdiag;
3423: return(0);
3424: }
3425: MatMarkDiagonal_SeqAIJ(A);
3426: if (!a->ibdiag) {
3427: PetscMalloc1(bs2*mbs,&a->ibdiag);
3428: PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));
3429: }
3430: diag = a->ibdiag;
3431: if (values) *values = a->ibdiag;
3432: /* factor and invert each block */
3433: switch (bs) {
3434: case 1:
3435: for (i=0; i<mbs; i++) {
3436: MatGetValues(A,1,&i,1,&i,diag+i);
3437: if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) {
3438: if (allowzeropivot) {
3439: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3440: A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]);
3441: A->factorerror_zeropivot_row = i;
3442: PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);
3443: } 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);
3444: }
3445: diag[i] = (PetscScalar)1.0 / (diag[i] + shift);
3446: }
3447: break;
3448: case 2:
3449: for (i=0; i<mbs; i++) {
3450: ij[0] = 2*i; ij[1] = 2*i + 1;
3451: MatGetValues(A,2,ij,2,ij,diag);
3452: PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
3453: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3454: PetscKernel_A_gets_transpose_A_2(diag);
3455: diag += 4;
3456: }
3457: break;
3458: case 3:
3459: for (i=0; i<mbs; i++) {
3460: ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2;
3461: MatGetValues(A,3,ij,3,ij,diag);
3462: PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
3463: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3464: PetscKernel_A_gets_transpose_A_3(diag);
3465: diag += 9;
3466: }
3467: break;
3468: case 4:
3469: for (i=0; i<mbs; i++) {
3470: ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3;
3471: MatGetValues(A,4,ij,4,ij,diag);
3472: PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
3473: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3474: PetscKernel_A_gets_transpose_A_4(diag);
3475: diag += 16;
3476: }
3477: break;
3478: case 5:
3479: for (i=0; i<mbs; i++) {
3480: ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4;
3481: MatGetValues(A,5,ij,5,ij,diag);
3482: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
3483: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3484: PetscKernel_A_gets_transpose_A_5(diag);
3485: diag += 25;
3486: }
3487: break;
3488: case 6:
3489: for (i=0; i<mbs; i++) {
3490: 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;
3491: MatGetValues(A,6,ij,6,ij,diag);
3492: PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
3493: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3494: PetscKernel_A_gets_transpose_A_6(diag);
3495: diag += 36;
3496: }
3497: break;
3498: case 7:
3499: for (i=0; i<mbs; i++) {
3500: 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;
3501: MatGetValues(A,7,ij,7,ij,diag);
3502: PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
3503: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3504: PetscKernel_A_gets_transpose_A_7(diag);
3505: diag += 49;
3506: }
3507: break;
3508: default:
3509: PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);
3510: for (i=0; i<mbs; i++) {
3511: for (j=0; j<bs; j++) {
3512: IJ[j] = bs*i + j;
3513: }
3514: MatGetValues(A,bs,IJ,bs,IJ,diag);
3515: PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
3516: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3517: PetscKernel_A_gets_transpose_A_N(diag,bs);
3518: diag += bs2;
3519: }
3520: PetscFree3(v_work,v_pivots,IJ);
3521: }
3522: a->ibdiagvalid = PETSC_TRUE;
3523: return(0);
3524: }
3526: static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx)
3527: {
3529: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3530: PetscScalar a;
3531: PetscInt m,n,i,j,col;
3534: if (!x->assembled) {
3535: MatGetSize(x,&m,&n);
3536: for (i=0; i<m; i++) {
3537: for (j=0; j<aij->imax[i]; j++) {
3538: PetscRandomGetValue(rctx,&a);
3539: col = (PetscInt)(n*PetscRealPart(a));
3540: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3541: }
3542: }
3543: } else {
3544: for (i=0; i<aij->nz; i++) {PetscRandomGetValue(rctx,aij->a+i);}
3545: }
3546: #if defined(PETSC_HAVE_DEVICE)
3547: if (x->offloadmask != PETSC_OFFLOAD_UNALLOCATED) x->offloadmask = PETSC_OFFLOAD_CPU;
3548: #endif
3549: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3550: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3551: return(0);
3552: }
3554: /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */
3555: PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx)
3556: {
3558: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3559: PetscScalar a;
3560: PetscInt m,n,i,j,col,nskip;
3563: nskip = high - low;
3564: MatGetSize(x,&m,&n);
3565: n -= nskip; /* shrink number of columns where nonzeros can be set */
3566: for (i=0; i<m; i++) {
3567: for (j=0; j<aij->imax[i]; j++) {
3568: PetscRandomGetValue(rctx,&a);
3569: col = (PetscInt)(n*PetscRealPart(a));
3570: if (col >= low) col += nskip; /* shift col rightward to skip the hole */
3571: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3572: }
3573: }
3574: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3575: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3576: return(0);
3577: }
3579: /* -------------------------------------------------------------------*/
3580: static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ,
3581: MatGetRow_SeqAIJ,
3582: MatRestoreRow_SeqAIJ,
3583: MatMult_SeqAIJ,
3584: /* 4*/ MatMultAdd_SeqAIJ,
3585: MatMultTranspose_SeqAIJ,
3586: MatMultTransposeAdd_SeqAIJ,
3587: NULL,
3588: NULL,
3589: NULL,
3590: /* 10*/ NULL,
3591: MatLUFactor_SeqAIJ,
3592: NULL,
3593: MatSOR_SeqAIJ,
3594: MatTranspose_SeqAIJ,
3595: /*1 5*/ MatGetInfo_SeqAIJ,
3596: MatEqual_SeqAIJ,
3597: MatGetDiagonal_SeqAIJ,
3598: MatDiagonalScale_SeqAIJ,
3599: MatNorm_SeqAIJ,
3600: /* 20*/ NULL,
3601: MatAssemblyEnd_SeqAIJ,
3602: MatSetOption_SeqAIJ,
3603: MatZeroEntries_SeqAIJ,
3604: /* 24*/ MatZeroRows_SeqAIJ,
3605: NULL,
3606: NULL,
3607: NULL,
3608: NULL,
3609: /* 29*/ MatSetUp_SeqAIJ,
3610: NULL,
3611: NULL,
3612: NULL,
3613: NULL,
3614: /* 34*/ MatDuplicate_SeqAIJ,
3615: NULL,
3616: NULL,
3617: MatILUFactor_SeqAIJ,
3618: NULL,
3619: /* 39*/ MatAXPY_SeqAIJ,
3620: MatCreateSubMatrices_SeqAIJ,
3621: MatIncreaseOverlap_SeqAIJ,
3622: MatGetValues_SeqAIJ,
3623: MatCopy_SeqAIJ,
3624: /* 44*/ MatGetRowMax_SeqAIJ,
3625: MatScale_SeqAIJ,
3626: MatShift_SeqAIJ,
3627: MatDiagonalSet_SeqAIJ,
3628: MatZeroRowsColumns_SeqAIJ,
3629: /* 49*/ MatSetRandom_SeqAIJ,
3630: MatGetRowIJ_SeqAIJ,
3631: MatRestoreRowIJ_SeqAIJ,
3632: MatGetColumnIJ_SeqAIJ,
3633: MatRestoreColumnIJ_SeqAIJ,
3634: /* 54*/ MatFDColoringCreate_SeqXAIJ,
3635: NULL,
3636: NULL,
3637: MatPermute_SeqAIJ,
3638: NULL,
3639: /* 59*/ NULL,
3640: MatDestroy_SeqAIJ,
3641: MatView_SeqAIJ,
3642: NULL,
3643: NULL,
3644: /* 64*/ NULL,
3645: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ,
3646: NULL,
3647: NULL,
3648: NULL,
3649: /* 69*/ MatGetRowMaxAbs_SeqAIJ,
3650: MatGetRowMinAbs_SeqAIJ,
3651: NULL,
3652: NULL,
3653: NULL,
3654: /* 74*/ NULL,
3655: MatFDColoringApply_AIJ,
3656: NULL,
3657: NULL,
3658: NULL,
3659: /* 79*/ MatFindZeroDiagonals_SeqAIJ,
3660: NULL,
3661: NULL,
3662: NULL,
3663: MatLoad_SeqAIJ,
3664: /* 84*/ MatIsSymmetric_SeqAIJ,
3665: MatIsHermitian_SeqAIJ,
3666: NULL,
3667: NULL,
3668: NULL,
3669: /* 89*/ NULL,
3670: NULL,
3671: MatMatMultNumeric_SeqAIJ_SeqAIJ,
3672: NULL,
3673: NULL,
3674: /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy,
3675: NULL,
3676: NULL,
3677: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ,
3678: NULL,
3679: /* 99*/ MatProductSetFromOptions_SeqAIJ,
3680: NULL,
3681: NULL,
3682: MatConjugate_SeqAIJ,
3683: NULL,
3684: /*104*/ MatSetValuesRow_SeqAIJ,
3685: MatRealPart_SeqAIJ,
3686: MatImaginaryPart_SeqAIJ,
3687: NULL,
3688: NULL,
3689: /*109*/ MatMatSolve_SeqAIJ,
3690: NULL,
3691: MatGetRowMin_SeqAIJ,
3692: NULL,
3693: MatMissingDiagonal_SeqAIJ,
3694: /*114*/ NULL,
3695: NULL,
3696: NULL,
3697: NULL,
3698: NULL,
3699: /*119*/ NULL,
3700: NULL,
3701: NULL,
3702: NULL,
3703: MatGetMultiProcBlock_SeqAIJ,
3704: /*124*/ MatFindNonzeroRows_SeqAIJ,
3705: MatGetColumnReductions_SeqAIJ,
3706: MatInvertBlockDiagonal_SeqAIJ,
3707: MatInvertVariableBlockDiagonal_SeqAIJ,
3708: NULL,
3709: /*129*/ NULL,
3710: NULL,
3711: NULL,
3712: MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ,
3713: MatTransposeColoringCreate_SeqAIJ,
3714: /*134*/ MatTransColoringApplySpToDen_SeqAIJ,
3715: MatTransColoringApplyDenToSp_SeqAIJ,
3716: NULL,
3717: NULL,
3718: MatRARtNumeric_SeqAIJ_SeqAIJ,
3719: /*139*/NULL,
3720: NULL,
3721: NULL,
3722: MatFDColoringSetUp_SeqXAIJ,
3723: MatFindOffBlockDiagonalEntries_SeqAIJ,
3724: MatCreateMPIMatConcatenateSeqMat_SeqAIJ,
3725: /*145*/MatDestroySubMatrices_SeqAIJ,
3726: NULL,
3727: NULL
3728: };
3730: PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices)
3731: {
3732: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3733: PetscInt i,nz,n;
3736: nz = aij->maxnz;
3737: n = mat->rmap->n;
3738: for (i=0; i<nz; i++) {
3739: aij->j[i] = indices[i];
3740: }
3741: aij->nz = nz;
3742: for (i=0; i<n; i++) {
3743: aij->ilen[i] = aij->imax[i];
3744: }
3745: return(0);
3746: }
3748: /*
3749: * Given a sparse matrix with global column indices, compact it by using a local column space.
3750: * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable()
3751: */
3752: PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping)
3753: {
3754: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3755: PetscTable gid1_lid1;
3756: PetscTablePosition tpos;
3757: PetscInt gid,lid,i,ec,nz = aij->nz;
3758: PetscInt *garray,*jj = aij->j;
3759: PetscErrorCode ierr;
3764: /* use a table */
3765: PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);
3766: ec = 0;
3767: for (i=0; i<nz; i++) {
3768: PetscInt data,gid1 = jj[i] + 1;
3769: PetscTableFind(gid1_lid1,gid1,&data);
3770: if (!data) {
3771: /* one based table */
3772: PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);
3773: }
3774: }
3775: /* form array of columns we need */
3776: PetscMalloc1(ec,&garray);
3777: PetscTableGetHeadPosition(gid1_lid1,&tpos);
3778: while (tpos) {
3779: PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);
3780: gid--;
3781: lid--;
3782: garray[lid] = gid;
3783: }
3784: PetscSortInt(ec,garray); /* sort, and rebuild */
3785: PetscTableRemoveAll(gid1_lid1);
3786: for (i=0; i<ec; i++) {
3787: PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);
3788: }
3789: /* compact out the extra columns in B */
3790: for (i=0; i<nz; i++) {
3791: PetscInt gid1 = jj[i] + 1;
3792: PetscTableFind(gid1_lid1,gid1,&lid);
3793: lid--;
3794: jj[i] = lid;
3795: }
3796: PetscLayoutDestroy(&mat->cmap);
3797: PetscTableDestroy(&gid1_lid1);
3798: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap);
3799: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);
3800: ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);
3801: return(0);
3802: }
3804: /*@
3805: MatSeqAIJSetColumnIndices - Set the column indices for all the rows
3806: in the matrix.
3808: Input Parameters:
3809: + mat - the SeqAIJ matrix
3810: - indices - the column indices
3812: Level: advanced
3814: Notes:
3815: This can be called if you have precomputed the nonzero structure of the
3816: matrix and want to provide it to the matrix object to improve the performance
3817: of the MatSetValues() operation.
3819: You MUST have set the correct numbers of nonzeros per row in the call to
3820: MatCreateSeqAIJ(), and the columns indices MUST be sorted.
3822: MUST be called before any calls to MatSetValues();
3824: The indices should start with zero, not one.
3826: @*/
3827: PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices)
3828: {
3834: PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
3835: return(0);
3836: }
3838: /* ----------------------------------------------------------------------------------------*/
3840: PetscErrorCode MatStoreValues_SeqAIJ(Mat mat)
3841: {
3842: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3844: size_t nz = aij->i[mat->rmap->n];
3847: if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3849: /* allocate space for values if not already there */
3850: if (!aij->saved_values) {
3851: PetscMalloc1(nz+1,&aij->saved_values);
3852: PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
3853: }
3855: /* copy values over */
3856: PetscArraycpy(aij->saved_values,aij->a,nz);
3857: return(0);
3858: }
3860: /*@
3861: MatStoreValues - Stashes a copy of the matrix values; this allows, for
3862: example, reuse of the linear part of a Jacobian, while recomputing the
3863: nonlinear portion.
3865: Collect on Mat
3867: Input Parameters:
3868: . mat - the matrix (currently only AIJ matrices support this option)
3870: Level: advanced
3872: Common Usage, with SNESSolve():
3873: $ Create Jacobian matrix
3874: $ Set linear terms into matrix
3875: $ Apply boundary conditions to matrix, at this time matrix must have
3876: $ final nonzero structure (i.e. setting the nonlinear terms and applying
3877: $ boundary conditions again will not change the nonzero structure
3878: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3879: $ MatStoreValues(mat);
3880: $ Call SNESSetJacobian() with matrix
3881: $ In your Jacobian routine
3882: $ MatRetrieveValues(mat);
3883: $ Set nonlinear terms in matrix
3885: Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself:
3886: $ // build linear portion of Jacobian
3887: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3888: $ MatStoreValues(mat);
3889: $ loop over nonlinear iterations
3890: $ MatRetrieveValues(mat);
3891: $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian
3892: $ // call MatAssemblyBegin/End() on matrix
3893: $ Solve linear system with Jacobian
3894: $ endloop
3896: Notes:
3897: Matrix must already be assemblied before calling this routine
3898: Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before
3899: calling this routine.
3901: When this is called multiple times it overwrites the previous set of stored values
3902: and does not allocated additional space.
3904: .seealso: MatRetrieveValues()
3906: @*/
3907: PetscErrorCode MatStoreValues(Mat mat)
3908: {
3913: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3914: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3915: PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));
3916: return(0);
3917: }
3919: PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat)
3920: {
3921: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3923: PetscInt nz = aij->i[mat->rmap->n];
3926: if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3927: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
3928: /* copy values over */
3929: PetscArraycpy(aij->a,aij->saved_values,nz);
3930: return(0);
3931: }
3933: /*@
3934: MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for
3935: example, reuse of the linear part of a Jacobian, while recomputing the
3936: nonlinear portion.
3938: Collect on Mat
3940: Input Parameters:
3941: . mat - the matrix (currently only AIJ matrices support this option)
3943: Level: advanced
3945: .seealso: MatStoreValues()
3947: @*/
3948: PetscErrorCode MatRetrieveValues(Mat mat)
3949: {
3954: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3955: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3956: PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));
3957: return(0);
3958: }
3960: /* --------------------------------------------------------------------------------*/
3961: /*@C
3962: MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format
3963: (the default parallel PETSc format). For good matrix assembly performance
3964: the user should preallocate the matrix storage by setting the parameter nz
3965: (or the array nnz). By setting these parameters accurately, performance
3966: during matrix assembly can be increased by more than a factor of 50.
3968: Collective
3970: Input Parameters:
3971: + comm - MPI communicator, set to PETSC_COMM_SELF
3972: . m - number of rows
3973: . n - number of columns
3974: . nz - number of nonzeros per row (same for all rows)
3975: - nnz - array containing the number of nonzeros in the various rows
3976: (possibly different for each row) or NULL
3978: Output Parameter:
3979: . A - the matrix
3981: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3982: MatXXXXSetPreallocation() paradigm instead of this routine directly.
3983: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
3985: Notes:
3986: If nnz is given then nz is ignored
3988: The AIJ format (also called the Yale sparse matrix format or
3989: compressed row storage), is fully compatible with standard Fortran 77
3990: storage. That is, the stored row and column indices can begin at
3991: either one (as in Fortran) or zero. See the users' manual for details.
3993: Specify the preallocated storage with either nz or nnz (not both).
3994: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3995: allocation. For large problems you MUST preallocate memory or you
3996: will get TERRIBLE performance, see the users' manual chapter on matrices.
3998: By default, this format uses inodes (identical nodes) when possible, to
3999: improve numerical efficiency of matrix-vector products and solves. We
4000: search for consecutive rows with the same nonzero structure, thereby
4001: reusing matrix information to achieve increased efficiency.
4003: Options Database Keys:
4004: + -mat_no_inode - Do not use inodes
4005: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
4007: Level: intermediate
4009: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays()
4011: @*/
4012: PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
4013: {
4017: MatCreate(comm,A);
4018: MatSetSizes(*A,m,n,m,n);
4019: MatSetType(*A,MATSEQAIJ);
4020: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
4021: return(0);
4022: }
4024: /*@C
4025: MatSeqAIJSetPreallocation - For good matrix assembly performance
4026: the user should preallocate the matrix storage by setting the parameter nz
4027: (or the array nnz). By setting these parameters accurately, performance
4028: during matrix assembly can be increased by more than a factor of 50.
4030: Collective
4032: Input Parameters:
4033: + B - The matrix
4034: . nz - number of nonzeros per row (same for all rows)
4035: - nnz - array containing the number of nonzeros in the various rows
4036: (possibly different for each row) or NULL
4038: Notes:
4039: If nnz is given then nz is ignored
4041: The AIJ format (also called the Yale sparse matrix format or
4042: compressed row storage), is fully compatible with standard Fortran 77
4043: storage. That is, the stored row and column indices can begin at
4044: either one (as in Fortran) or zero. See the users' manual for details.
4046: Specify the preallocated storage with either nz or nnz (not both).
4047: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
4048: allocation. For large problems you MUST preallocate memory or you
4049: will get TERRIBLE performance, see the users' manual chapter on matrices.
4051: You can call MatGetInfo() to get information on how effective the preallocation was;
4052: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
4053: You can also run with the option -info and look for messages with the string
4054: malloc in them to see if additional memory allocation was needed.
4056: Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix
4057: entries or columns indices
4059: By default, this format uses inodes (identical nodes) when possible, to
4060: improve numerical efficiency of matrix-vector products and solves. We
4061: search for consecutive rows with the same nonzero structure, thereby
4062: reusing matrix information to achieve increased efficiency.
4064: Options Database Keys:
4065: + -mat_no_inode - Do not use inodes
4066: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
4068: Level: intermediate
4070: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo(),
4071: MatSeqAIJSetTotalPreallocation()
4073: @*/
4074: PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[])
4075: {
4081: PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));
4082: return(0);
4083: }
4085: PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz)
4086: {
4087: Mat_SeqAIJ *b;
4088: PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;
4090: PetscInt i;
4093: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
4094: if (nz == MAT_SKIP_ALLOCATION) {
4095: skipallocation = PETSC_TRUE;
4096: nz = 0;
4097: }
4098: PetscLayoutSetUp(B->rmap);
4099: PetscLayoutSetUp(B->cmap);
4101: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
4102: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
4103: if (PetscUnlikelyDebug(nnz)) {
4104: for (i=0; i<B->rmap->n; i++) {
4105: 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]);
4106: 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);
4107: }
4108: }
4110: B->preallocated = PETSC_TRUE;
4112: b = (Mat_SeqAIJ*)B->data;
4114: if (!skipallocation) {
4115: if (!b->imax) {
4116: PetscMalloc1(B->rmap->n,&b->imax);
4117: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
4118: }
4119: if (!b->ilen) {
4120: /* b->ilen will count nonzeros in each row so far. */
4121: PetscCalloc1(B->rmap->n,&b->ilen);
4122: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
4123: } else {
4124: PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt));
4125: }
4126: if (!b->ipre) {
4127: PetscMalloc1(B->rmap->n,&b->ipre);
4128: PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));
4129: }
4130: if (!nnz) {
4131: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10;
4132: else if (nz < 0) nz = 1;
4133: nz = PetscMin(nz,B->cmap->n);
4134: for (i=0; i<B->rmap->n; i++) b->imax[i] = nz;
4135: nz = nz*B->rmap->n;
4136: } else {
4137: PetscInt64 nz64 = 0;
4138: for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
4139: PetscIntCast(nz64,&nz);
4140: }
4142: /* allocate the matrix space */
4143: /* FIXME: should B's old memory be unlogged? */
4144: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
4145: if (B->structure_only) {
4146: PetscMalloc1(nz,&b->j);
4147: PetscMalloc1(B->rmap->n+1,&b->i);
4148: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));
4149: } else {
4150: PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);
4151: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
4152: }
4153: b->i[0] = 0;
4154: for (i=1; i<B->rmap->n+1; i++) {
4155: b->i[i] = b->i[i-1] + b->imax[i-1];
4156: }
4157: if (B->structure_only) {
4158: b->singlemalloc = PETSC_FALSE;
4159: b->free_a = PETSC_FALSE;
4160: } else {
4161: b->singlemalloc = PETSC_TRUE;
4162: b->free_a = PETSC_TRUE;
4163: }
4164: b->free_ij = PETSC_TRUE;
4165: } else {
4166: b->free_a = PETSC_FALSE;
4167: b->free_ij = PETSC_FALSE;
4168: }
4170: if (b->ipre && nnz != b->ipre && b->imax) {
4171: /* reserve user-requested sparsity */
4172: PetscArraycpy(b->ipre,b->imax,B->rmap->n);
4173: }
4175: b->nz = 0;
4176: b->maxnz = nz;
4177: B->info.nz_unneeded = (double)b->maxnz;
4178: if (realalloc) {
4179: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
4180: }
4181: B->was_assembled = PETSC_FALSE;
4182: B->assembled = PETSC_FALSE;
4183: return(0);
4184: }
4186: PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A)
4187: {
4188: Mat_SeqAIJ *a;
4189: PetscInt i;
4195: /* Check local size. If zero, then return */
4196: if (!A->rmap->n) return(0);
4198: a = (Mat_SeqAIJ*)A->data;
4199: /* if no saved info, we error out */
4200: if (!a->ipre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info \n");
4202: 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");
4204: PetscArraycpy(a->imax,a->ipre,A->rmap->n);
4205: PetscArrayzero(a->ilen,A->rmap->n);
4206: a->i[0] = 0;
4207: for (i=1; i<A->rmap->n+1; i++) {
4208: a->i[i] = a->i[i-1] + a->imax[i-1];
4209: }
4210: A->preallocated = PETSC_TRUE;
4211: a->nz = 0;
4212: a->maxnz = a->i[A->rmap->n];
4213: A->info.nz_unneeded = (double)a->maxnz;
4214: A->was_assembled = PETSC_FALSE;
4215: A->assembled = PETSC_FALSE;
4216: return(0);
4217: }
4219: /*@
4220: MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format.
4222: Input Parameters:
4223: + B - the matrix
4224: . i - the indices into j for the start of each row (starts with zero)
4225: . j - the column indices for each row (starts with zero) these must be sorted for each row
4226: - v - optional values in the matrix
4228: Level: developer
4230: Notes:
4231: The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays()
4233: This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero
4234: structure will be the union of all the previous nonzero structures.
4236: Developer Notes:
4237: 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
4238: then just copies the v values directly with PetscMemcpy().
4240: This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them.
4242: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation()
4243: @*/
4244: PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[])
4245: {
4251: PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));
4252: return(0);
4253: }
4255: PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[])
4256: {
4257: PetscInt i;
4258: PetscInt m,n;
4259: PetscInt nz;
4260: PetscInt *nnz;
4264: if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]);
4266: PetscLayoutSetUp(B->rmap);
4267: PetscLayoutSetUp(B->cmap);
4269: MatGetSize(B, &m, &n);
4270: PetscMalloc1(m+1, &nnz);
4271: for (i = 0; i < m; i++) {
4272: nz = Ii[i+1]- Ii[i];
4273: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz);
4274: nnz[i] = nz;
4275: }
4276: MatSeqAIJSetPreallocation(B, 0, nnz);
4277: PetscFree(nnz);
4279: for (i = 0; i < m; i++) {
4280: MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES);
4281: }
4283: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
4284: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
4286: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
4287: return(0);
4288: }
4290: /*@
4291: MatSeqAIJKron - Computes C, the Kronecker product of A and B.
4293: Input Parameters:
4294: + A - left-hand side matrix
4295: . B - right-hand side matrix
4296: - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4298: Output Parameter:
4299: . C - Kronecker product of A and B
4301: Level: intermediate
4303: Notes:
4304: MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron().
4306: .seealso: MatCreateSeqAIJ(), MATSEQAIJ, MATKAIJ, MatReuse
4307: @*/
4308: PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C)
4309: {
4318: if (reuse == MAT_REUSE_MATRIX) {
4321: }
4322: PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C));
4323: return(0);
4324: }
4326: PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C)
4327: {
4328: Mat newmat;
4329: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4330: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
4331: PetscScalar *v;
4332: 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;
4333: PetscBool flg;
4337: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4338: if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4339: if (B->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4340: if (!B->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4341: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg);
4342: if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name);
4343: if (reuse != MAT_INITIAL_MATRIX && reuse != MAT_REUSE_MATRIX) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse);
4344: if (reuse == MAT_INITIAL_MATRIX) {
4345: PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j);
4346: MatCreate(PETSC_COMM_SELF,&newmat);
4347: MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn);
4348: MatSetType(newmat,MATAIJ);
4349: i[0] = 0;
4350: for (m = 0; m < am; ++m) {
4351: for (p = 0; p < bm; ++p) {
4352: i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]);
4353: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4354: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4355: j[nnz++] = a->j[n]*bn + b->j[q];
4356: }
4357: }
4358: }
4359: }
4360: MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL);
4361: *C = newmat;
4362: PetscFree2(i,j);
4363: nnz = 0;
4364: }
4365: MatSeqAIJGetArray(*C,&v);
4366: for (m = 0; m < am; ++m) {
4367: for (p = 0; p < bm; ++p) {
4368: for (n = a->i[m]; n < a->i[m+1]; ++n) {
4369: for (q = b->i[p]; q < b->i[p+1]; ++q) {
4370: v[nnz++] = a->a[n] * b->a[q];
4371: }
4372: }
4373: }
4374: }
4375: MatSeqAIJRestoreArray(*C,&v);
4376: return(0);
4377: }
4379: #include <../src/mat/impls/dense/seq/dense.h>
4380: #include <petsc/private/kernels/petscaxpy.h>
4382: /*
4383: Computes (B'*A')' since computing B*A directly is untenable
4385: n p p
4386: [ ] [ ] [ ]
4387: m [ A ] * n [ B ] = m [ C ]
4388: [ ] [ ] [ ]
4390: */
4391: PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C)
4392: {
4393: PetscErrorCode ierr;
4394: Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data;
4395: Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data;
4396: Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data;
4397: PetscInt i,j,n,m,q,p;
4398: const PetscInt *ii,*idx;
4399: const PetscScalar *b,*a,*a_q;
4400: PetscScalar *c,*c_q;
4401: PetscInt clda = sub_c->lda;
4402: PetscInt alda = sub_a->lda;
4405: m = A->rmap->n;
4406: n = A->cmap->n;
4407: p = B->cmap->n;
4408: a = sub_a->v;
4409: b = sub_b->a;
4410: c = sub_c->v;
4411: if (clda == m) {
4412: PetscArrayzero(c,m*p);
4413: } else {
4414: for (j=0;j<p;j++)
4415: for (i=0;i<m;i++)
4416: c[j*clda + i] = 0.0;
4417: }
4418: ii = sub_b->i;
4419: idx = sub_b->j;
4420: for (i=0; i<n; i++) {
4421: q = ii[i+1] - ii[i];
4422: while (q-->0) {
4423: c_q = c + clda*(*idx);
4424: a_q = a + alda*i;
4425: PetscKernelAXPY(c_q,*b,a_q,m);
4426: idx++;
4427: b++;
4428: }
4429: }
4430: return(0);
4431: }
4433: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C)
4434: {
4436: PetscInt m=A->rmap->n,n=B->cmap->n;
4437: PetscBool cisdense;
4440: 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);
4441: MatSetSizes(C,m,n,m,n);
4442: MatSetBlockSizesFromMats(C,A,B);
4443: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
4444: if (!cisdense) {
4445: MatSetType(C,MATDENSE);
4446: }
4447: MatSetUp(C);
4449: C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ;
4450: return(0);
4451: }
4453: /* ----------------------------------------------------------------*/
4454: /*MC
4455: MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices,
4456: based on compressed sparse row format.
4458: Options Database Keys:
4459: . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions()
4461: Level: beginner
4463: Notes:
4464: MatSetValues() may be called for this matrix type with a NULL argument for the numerical values,
4465: in this case the values associated with the rows and columns one passes in are set to zero
4466: in the matrix
4468: MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no
4469: space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored
4471: Developer Notes:
4472: It would be nice if all matrix formats supported passing NULL in for the numerical values
4474: .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType
4475: M*/
4477: /*MC
4478: MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices.
4480: This matrix type is identical to MATSEQAIJ when constructed with a single process communicator,
4481: and MATMPIAIJ otherwise. As a result, for single process communicators,
4482: MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation() is supported
4483: for communicators controlling multiple processes. It is recommended that you call both of
4484: the above preallocation routines for simplicity.
4486: Options Database Keys:
4487: . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions()
4489: Developer Notes:
4490: Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when
4491: enough exist.
4493: Level: beginner
4495: .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ
4496: M*/
4498: /*MC
4499: MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices.
4501: This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator,
4502: and MATMPIAIJCRL otherwise. As a result, for single process communicators,
4503: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
4504: for communicators controlling multiple processes. It is recommended that you call both of
4505: the above preallocation routines for simplicity.
4507: Options Database Keys:
4508: . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions()
4510: Level: beginner
4512: .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL
4513: M*/
4515: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*);
4516: #if defined(PETSC_HAVE_ELEMENTAL)
4517: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
4518: #endif
4519: #if defined(PETSC_HAVE_SCALAPACK)
4520: PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
4521: #endif
4522: #if defined(PETSC_HAVE_HYPRE)
4523: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*);
4524: #endif
4526: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*);
4527: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*);
4528: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat);
4530: /*@C
4531: MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored
4533: Not Collective
4535: Input Parameter:
4536: . mat - a MATSEQAIJ matrix
4538: Output Parameter:
4539: . array - pointer to the data
4541: Level: intermediate
4543: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4544: @*/
4545: PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array)
4546: {
4550: PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));
4551: #if defined(PETSC_HAVE_DEVICE)
4552: if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED) A->offloadmask = PETSC_OFFLOAD_CPU;
4553: #endif
4554: return(0);
4555: }
4557: /*@C
4558: MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored
4560: Not Collective
4562: Input Parameter:
4563: . mat - a MATSEQAIJ matrix
4565: Output Parameter:
4566: . array - pointer to the data
4568: Level: intermediate
4570: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead()
4571: @*/
4572: PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array)
4573: {
4574: #if defined(PETSC_HAVE_DEVICE)
4575: PetscOffloadMask oval;
4576: #endif
4580: #if defined(PETSC_HAVE_DEVICE)
4581: oval = A->offloadmask;
4582: #endif
4583: MatSeqAIJGetArray(A,(PetscScalar**)array);
4584: #if defined(PETSC_HAVE_DEVICE)
4585: if (oval == PETSC_OFFLOAD_GPU || oval == PETSC_OFFLOAD_BOTH) A->offloadmask = PETSC_OFFLOAD_BOTH;
4586: #endif
4587: return(0);
4588: }
4590: /*@C
4591: MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead
4593: Not Collective
4595: Input Parameter:
4596: . mat - a MATSEQAIJ matrix
4598: Output Parameter:
4599: . array - pointer to the data
4601: Level: intermediate
4603: .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead()
4604: @*/
4605: PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array)
4606: {
4607: #if defined(PETSC_HAVE_DEVICE)
4608: PetscOffloadMask oval;
4609: #endif
4613: #if defined(PETSC_HAVE_DEVICE)
4614: oval = A->offloadmask;
4615: #endif
4616: MatSeqAIJRestoreArray(A,(PetscScalar**)array);
4617: #if defined(PETSC_HAVE_DEVICE)
4618: A->offloadmask = oval;
4619: #endif
4620: return(0);
4621: }
4623: /*@C
4624: MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row
4626: Not Collective
4628: Input Parameter:
4629: . mat - a MATSEQAIJ matrix
4631: Output Parameter:
4632: . nz - the maximum number of nonzeros in any row
4634: Level: intermediate
4636: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
4637: @*/
4638: PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz)
4639: {
4640: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
4643: *nz = aij->rmax;
4644: return(0);
4645: }
4647: /*@C
4648: MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray()
4650: Not Collective
4652: Input Parameters:
4653: + mat - a MATSEQAIJ matrix
4654: - array - pointer to the data
4656: Level: intermediate
4658: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90()
4659: @*/
4660: PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array)
4661: {
4665: PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
4666: return(0);
4667: }
4669: #if defined(PETSC_HAVE_CUDA)
4670: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*);
4671: #endif
4672: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4673: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*);
4674: #endif
4676: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B)
4677: {
4678: Mat_SeqAIJ *b;
4680: PetscMPIInt size;
4683: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
4684: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1");
4686: PetscNewLog(B,&b);
4688: B->data = (void*)b;
4690: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
4691: if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull;
4693: b->row = NULL;
4694: b->col = NULL;
4695: b->icol = NULL;
4696: b->reallocs = 0;
4697: b->ignorezeroentries = PETSC_FALSE;
4698: b->roworiented = PETSC_TRUE;
4699: b->nonew = 0;
4700: b->diag = NULL;
4701: b->solve_work = NULL;
4702: B->spptr = NULL;
4703: b->saved_values = NULL;
4704: b->idiag = NULL;
4705: b->mdiag = NULL;
4706: b->ssor_work = NULL;
4707: b->omega = 1.0;
4708: b->fshift = 0.0;
4709: b->idiagvalid = PETSC_FALSE;
4710: b->ibdiagvalid = PETSC_FALSE;
4711: b->keepnonzeropattern = PETSC_FALSE;
4713: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4714: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);
4715: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);
4717: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4718: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);
4719: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);
4720: #endif
4722: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);
4723: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);
4724: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);
4725: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);
4726: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);
4727: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);
4728: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);
4729: #if defined(PETSC_HAVE_MKL_SPARSE)
4730: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);
4731: #endif
4732: #if defined(PETSC_HAVE_CUDA)
4733: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);
4734: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4735: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ);
4736: #endif
4737: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
4738: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos);
4739: #endif
4740: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);
4741: #if defined(PETSC_HAVE_ELEMENTAL)
4742: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);
4743: #endif
4744: #if defined(PETSC_HAVE_SCALAPACK)
4745: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK);
4746: #endif
4747: #if defined(PETSC_HAVE_HYPRE)
4748: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);
4749: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ);
4750: #endif
4751: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);
4752: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);
4753: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);
4754: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);
4755: PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);
4756: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);
4757: PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);
4758: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);
4759: PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);
4760: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ);
4761: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ);
4762: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ);
4763: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ);
4764: MatCreate_SeqAIJ_Inode(B);
4765: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
4766: MatSeqAIJSetTypeFromOptions(B); /* this allows changing the matrix subtype to say MATSEQAIJPERM */
4767: return(0);
4768: }
4770: /*
4771: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
4772: */
4773: PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
4774: {
4775: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data;
4777: PetscInt m = A->rmap->n,i;
4780: if (!A->assembled && cpvalues!=MAT_DO_NOT_COPY_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix");
4782: C->factortype = A->factortype;
4783: c->row = NULL;
4784: c->col = NULL;
4785: c->icol = NULL;
4786: c->reallocs = 0;
4788: C->assembled = PETSC_TRUE;
4790: PetscLayoutReference(A->rmap,&C->rmap);
4791: PetscLayoutReference(A->cmap,&C->cmap);
4793: PetscMalloc1(m,&c->imax);
4794: PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt));
4795: PetscMalloc1(m,&c->ilen);
4796: PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt));
4797: PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));
4799: /* allocate the matrix space */
4800: if (mallocmatspace) {
4801: PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);
4802: PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));
4804: c->singlemalloc = PETSC_TRUE;
4806: PetscArraycpy(c->i,a->i,m+1);
4807: if (m > 0) {
4808: PetscArraycpy(c->j,a->j,a->i[m]);
4809: if (cpvalues == MAT_COPY_VALUES) {
4810: const PetscScalar *aa;
4812: MatSeqAIJGetArrayRead(A,&aa);
4813: PetscArraycpy(c->a,aa,a->i[m]);
4814: MatSeqAIJGetArrayRead(A,&aa);
4815: } else {
4816: PetscArrayzero(c->a,a->i[m]);
4817: }
4818: }
4819: }
4821: c->ignorezeroentries = a->ignorezeroentries;
4822: c->roworiented = a->roworiented;
4823: c->nonew = a->nonew;
4824: if (a->diag) {
4825: PetscMalloc1(m+1,&c->diag);
4826: PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt));
4827: PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));
4828: } else c->diag = NULL;
4830: c->solve_work = NULL;
4831: c->saved_values = NULL;
4832: c->idiag = NULL;
4833: c->ssor_work = NULL;
4834: c->keepnonzeropattern = a->keepnonzeropattern;
4835: c->free_a = PETSC_TRUE;
4836: c->free_ij = PETSC_TRUE;
4838: c->rmax = a->rmax;
4839: c->nz = a->nz;
4840: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
4841: C->preallocated = PETSC_TRUE;
4843: c->compressedrow.use = a->compressedrow.use;
4844: c->compressedrow.nrows = a->compressedrow.nrows;
4845: if (a->compressedrow.use) {
4846: i = a->compressedrow.nrows;
4847: PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);
4848: PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);
4849: PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);
4850: } else {
4851: c->compressedrow.use = PETSC_FALSE;
4852: c->compressedrow.i = NULL;
4853: c->compressedrow.rindex = NULL;
4854: }
4855: c->nonzerorowcnt = a->nonzerorowcnt;
4856: C->nonzerostate = A->nonzerostate;
4858: MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);
4859: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
4860: return(0);
4861: }
4863: PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
4864: {
4868: MatCreate(PetscObjectComm((PetscObject)A),B);
4869: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
4870: if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) {
4871: MatSetBlockSizesFromMats(*B,A,A);
4872: }
4873: MatSetType(*B,((PetscObject)A)->type_name);
4874: MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);
4875: return(0);
4876: }
4878: PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer)
4879: {
4880: PetscBool isbinary, ishdf5;
4886: /* force binary viewer to load .info file if it has not yet done so */
4887: PetscViewerSetUp(viewer);
4888: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
4889: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
4890: if (isbinary) {
4891: MatLoad_SeqAIJ_Binary(newMat,viewer);
4892: } else if (ishdf5) {
4893: #if defined(PETSC_HAVE_HDF5)
4894: MatLoad_AIJ_HDF5(newMat,viewer);
4895: #else
4896: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
4897: #endif
4898: } else {
4899: 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);
4900: }
4901: return(0);
4902: }
4904: PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer)
4905: {
4906: Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data;
4908: PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i;
4911: PetscViewerSetUp(viewer);
4913: /* read in matrix header */
4914: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
4915: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
4916: M = header[1]; N = header[2]; nz = header[3];
4917: if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
4918: if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
4919: if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ");
4921: /* set block sizes from the viewer's .info file */
4922: MatLoad_Binary_BlockSizes(mat,viewer);
4923: /* set local and global sizes if not set already */
4924: if (mat->rmap->n < 0) mat->rmap->n = M;
4925: if (mat->cmap->n < 0) mat->cmap->n = N;
4926: if (mat->rmap->N < 0) mat->rmap->N = M;
4927: if (mat->cmap->N < 0) mat->cmap->N = N;
4928: PetscLayoutSetUp(mat->rmap);
4929: PetscLayoutSetUp(mat->cmap);
4931: /* check if the matrix sizes are correct */
4932: MatGetSize(mat,&rows,&cols);
4933: 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);
4935: /* read in row lengths */
4936: PetscMalloc1(M,&rowlens);
4937: PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT);
4938: /* check if sum(rowlens) is same as nz */
4939: sum = 0; for (i=0; i<M; i++) sum += rowlens[i];
4940: 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);
4941: /* preallocate and check sizes */
4942: MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens);
4943: MatGetSize(mat,&rows,&cols);
4944: 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);
4945: /* store row lengths */
4946: PetscArraycpy(a->ilen,rowlens,M);
4947: PetscFree(rowlens);
4949: /* fill in "i" row pointers */
4950: a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i];
4951: /* read in "j" column indices */
4952: PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT);
4953: /* read in "a" nonzero values */
4954: PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR);
4956: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
4957: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
4958: return(0);
4959: }
4961: PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg)
4962: {
4963: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data;
4965: #if defined(PETSC_USE_COMPLEX)
4966: PetscInt k;
4967: #endif
4970: /* If the matrix dimensions are not equal,or no of nonzeros */
4971: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) {
4972: *flg = PETSC_FALSE;
4973: return(0);
4974: }
4976: /* if the a->i are the same */
4977: PetscArraycmp(a->i,b->i,A->rmap->n+1,flg);
4978: if (!*flg) return(0);
4980: /* if a->j are the same */
4981: PetscArraycmp(a->j,b->j,a->nz,flg);
4982: if (!*flg) return(0);
4984: /* if a->a are the same */
4985: #if defined(PETSC_USE_COMPLEX)
4986: for (k=0; k<a->nz; k++) {
4987: if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) {
4988: *flg = PETSC_FALSE;
4989: return(0);
4990: }
4991: }
4992: #else
4993: PetscArraycmp(a->a,b->a,a->nz,flg);
4994: #endif
4995: return(0);
4996: }
4998: /*@
4999: MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format)
5000: provided by the user.
5002: Collective
5004: Input Parameters:
5005: + comm - must be an MPI communicator of size 1
5006: . m - number of rows
5007: . n - number of columns
5008: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix
5009: . j - column indices
5010: - a - matrix values
5012: Output Parameter:
5013: . mat - the matrix
5015: Level: intermediate
5017: Notes:
5018: The i, j, and a arrays are not copied by this routine, the user must free these arrays
5019: once the matrix is destroyed and not before
5021: You cannot set new nonzero locations into this matrix, that will generate an error.
5023: The i and j indices are 0 based
5025: The format which is used for the sparse matrix input, is equivalent to a
5026: row-major ordering.. i.e for the following matrix, the input data expected is
5027: as shown
5029: $ 1 0 0
5030: $ 2 0 3
5031: $ 4 5 6
5032: $
5033: $ i = {0,1,3,6} [size = nrow+1 = 3+1]
5034: $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row
5035: $ v = {1,2,3,4,5,6} [size = 6]
5037: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5039: @*/
5040: PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
5041: {
5043: PetscInt ii;
5044: Mat_SeqAIJ *aij;
5045: PetscInt jj;
5048: if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
5049: MatCreate(comm,mat);
5050: MatSetSizes(*mat,m,n,m,n);
5051: /* MatSetBlockSizes(*mat,,); */
5052: MatSetType(*mat,MATSEQAIJ);
5053: MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL);
5054: aij = (Mat_SeqAIJ*)(*mat)->data;
5055: PetscMalloc1(m,&aij->imax);
5056: PetscMalloc1(m,&aij->ilen);
5058: aij->i = i;
5059: aij->j = j;
5060: aij->a = a;
5061: aij->singlemalloc = PETSC_FALSE;
5062: aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
5063: aij->free_a = PETSC_FALSE;
5064: aij->free_ij = PETSC_FALSE;
5066: for (ii=0; ii<m; ii++) {
5067: aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
5068: if (PetscDefined(USE_DEBUG)) {
5069: 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]);
5070: for (jj=i[ii]+1; jj<i[ii+1]; jj++) {
5071: 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);
5072: 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);
5073: }
5074: }
5075: }
5076: if (PetscDefined(USE_DEBUG)) {
5077: for (ii=0; ii<aij->i[m]; ii++) {
5078: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]);
5079: 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]);
5080: }
5081: }
5083: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5084: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5085: return(0);
5086: }
5087: /*@C
5088: MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format)
5089: provided by the user.
5091: Collective
5093: Input Parameters:
5094: + comm - must be an MPI communicator of size 1
5095: . m - number of rows
5096: . n - number of columns
5097: . i - row indices
5098: . j - column indices
5099: . a - matrix values
5100: . nz - number of nonzeros
5101: - idx - 0 or 1 based
5103: Output Parameter:
5104: . mat - the matrix
5106: Level: intermediate
5108: Notes:
5109: 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,
5110: the input data expected is as shown
5111: .vb
5112: 1 0 0
5113: 2 0 3
5114: 4 5 6
5116: i = {0,1,1,2,2,2}
5117: j = {0,0,2,0,1,2}
5118: v = {1,2,3,4,5,6}
5119: .ve
5121: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
5123: @*/
5124: PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx)
5125: {
5127: PetscInt ii, *nnz, one = 1,row,col;
5130: PetscCalloc1(m,&nnz);
5131: for (ii = 0; ii < nz; ii++) {
5132: nnz[i[ii] - !!idx] += 1;
5133: }
5134: MatCreate(comm,mat);
5135: MatSetSizes(*mat,m,n,m,n);
5136: MatSetType(*mat,MATSEQAIJ);
5137: MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);
5138: for (ii = 0; ii < nz; ii++) {
5139: if (idx) {
5140: row = i[ii] - 1;
5141: col = j[ii] - 1;
5142: } else {
5143: row = i[ii];
5144: col = j[ii];
5145: }
5146: MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);
5147: }
5148: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
5149: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
5150: PetscFree(nnz);
5151: return(0);
5152: }
5154: PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A)
5155: {
5156: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
5160: a->idiagvalid = PETSC_FALSE;
5161: a->ibdiagvalid = PETSC_FALSE;
5163: MatSeqAIJInvalidateDiagonal_Inode(A);
5164: return(0);
5165: }
5167: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
5168: {
5170: PetscMPIInt size;
5173: MPI_Comm_size(comm,&size);
5174: if (size == 1) {
5175: if (scall == MAT_INITIAL_MATRIX) {
5176: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
5177: } else {
5178: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
5179: }
5180: } else {
5181: MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);
5182: }
5183: return(0);
5184: }
5186: /*
5187: Permute A into C's *local* index space using rowemb,colemb.
5188: The embedding are supposed to be injections and the above implies that the range of rowemb is a subset
5189: of [0,m), colemb is in [0,n).
5190: If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A.
5191: */
5192: PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B)
5193: {
5194: /* If making this function public, change the error returned in this function away from _PLIB. */
5196: Mat_SeqAIJ *Baij;
5197: PetscBool seqaij;
5198: PetscInt m,n,*nz,i,j,count;
5199: PetscScalar v;
5200: const PetscInt *rowindices,*colindices;
5203: if (!B) return(0);
5204: /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */
5205: PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
5206: if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type");
5207: if (rowemb) {
5208: ISGetLocalSize(rowemb,&m);
5209: 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);
5210: } else {
5211: if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix");
5212: }
5213: if (colemb) {
5214: ISGetLocalSize(colemb,&n);
5215: 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);
5216: } else {
5217: if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix");
5218: }
5220: Baij = (Mat_SeqAIJ*)(B->data);
5221: if (pattern == DIFFERENT_NONZERO_PATTERN) {
5222: PetscMalloc1(B->rmap->n,&nz);
5223: for (i=0; i<B->rmap->n; i++) {
5224: nz[i] = Baij->i[i+1] - Baij->i[i];
5225: }
5226: MatSeqAIJSetPreallocation(C,0,nz);
5227: PetscFree(nz);
5228: }
5229: if (pattern == SUBSET_NONZERO_PATTERN) {
5230: MatZeroEntries(C);
5231: }
5232: count = 0;
5233: rowindices = NULL;
5234: colindices = NULL;
5235: if (rowemb) {
5236: ISGetIndices(rowemb,&rowindices);
5237: }
5238: if (colemb) {
5239: ISGetIndices(colemb,&colindices);
5240: }
5241: for (i=0; i<B->rmap->n; i++) {
5242: PetscInt row;
5243: row = i;
5244: if (rowindices) row = rowindices[i];
5245: for (j=Baij->i[i]; j<Baij->i[i+1]; j++) {
5246: PetscInt col;
5247: col = Baij->j[count];
5248: if (colindices) col = colindices[col];
5249: v = Baij->a[count];
5250: MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);
5251: ++count;
5252: }
5253: }
5254: /* FIXME: set C's nonzerostate correctly. */
5255: /* Assembly for C is necessary. */
5256: C->preallocated = PETSC_TRUE;
5257: C->assembled = PETSC_TRUE;
5258: C->was_assembled = PETSC_FALSE;
5259: return(0);
5260: }
5262: PetscFunctionList MatSeqAIJList = NULL;
5264: /*@C
5265: MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype
5267: Collective on Mat
5269: Input Parameters:
5270: + mat - the matrix object
5271: - matype - matrix type
5273: Options Database Key:
5274: . -mat_seqai_type <method> - for example seqaijcrl
5276: Level: intermediate
5278: .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat
5279: @*/
5280: PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype)
5281: {
5282: PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*);
5283: PetscBool sametype;
5287: PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);
5288: if (sametype) return(0);
5290: PetscFunctionListFind(MatSeqAIJList,matype,&r);
5291: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype);
5292: (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);
5293: return(0);
5294: }
5296: /*@C
5297: MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices
5299: Not Collective
5301: Input Parameters:
5302: + name - name of a new user-defined matrix type, for example MATSEQAIJCRL
5303: - function - routine to convert to subtype
5305: Notes:
5306: MatSeqAIJRegister() may be called multiple times to add several user-defined solvers.
5308: Then, your matrix can be chosen with the procedural interface at runtime via the option
5309: $ -mat_seqaij_type my_mat
5311: Level: advanced
5313: .seealso: MatSeqAIJRegisterAll()
5315: Level: advanced
5316: @*/
5317: PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *))
5318: {
5322: MatInitializePackage();
5323: PetscFunctionListAdd(&MatSeqAIJList,sname,function);
5324: return(0);
5325: }
5327: PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE;
5329: /*@C
5330: MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ
5332: Not Collective
5334: Level: advanced
5336: .seealso: MatRegisterAll(), MatSeqAIJRegister()
5337: @*/
5338: PetscErrorCode MatSeqAIJRegisterAll(void)
5339: {
5343: if (MatSeqAIJRegisterAllCalled) return(0);
5344: MatSeqAIJRegisterAllCalled = PETSC_TRUE;
5346: MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);
5347: MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);
5348: MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);
5349: #if defined(PETSC_HAVE_MKL_SPARSE)
5350: MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);
5351: #endif
5352: #if defined(PETSC_HAVE_CUDA)
5353: MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE);
5354: #endif
5355: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
5356: MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos);
5357: #endif
5358: #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA)
5359: MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);
5360: #endif
5361: return(0);
5362: }
5364: /*
5365: Special version for direct calls from Fortran
5366: */
5367: #include <petsc/private/fortranimpl.h>
5368: #if defined(PETSC_HAVE_FORTRAN_CAPS)
5369: #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ
5370: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
5371: #define matsetvaluesseqaij_ matsetvaluesseqaij
5372: #endif
5374: /* Change these macros so can be used in void function */
5375: #undef CHKERRQ
5376: #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr)
5377: #undef SETERRQ2
5378: #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr)
5379: #undef SETERRQ3
5380: #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr)
5382: PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr)
5383: {
5384: Mat A = *AA;
5385: PetscInt m = *mm, n = *nn;
5386: InsertMode is = *isis;
5387: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
5388: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
5389: PetscInt *imax,*ai,*ailen;
5391: PetscInt *aj,nonew = a->nonew,lastcol = -1;
5392: MatScalar *ap,value,*aa;
5393: PetscBool ignorezeroentries = a->ignorezeroentries;
5394: PetscBool roworiented = a->roworiented;
5397: MatCheckPreallocated(A,1);
5398: imax = a->imax;
5399: ai = a->i;
5400: ailen = a->ilen;
5401: aj = a->j;
5402: aa = a->a;
5404: for (k=0; k<m; k++) { /* loop over added rows */
5405: row = im[k];
5406: if (row < 0) continue;
5407: if (PetscUnlikelyDebug(row >= A->rmap->n)) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
5408: rp = aj + ai[row]; ap = aa + ai[row];
5409: rmax = imax[row]; nrow = ailen[row];
5410: low = 0;
5411: high = nrow;
5412: for (l=0; l<n; l++) { /* loop over added columns */
5413: if (in[l] < 0) continue;
5414: if (PetscUnlikelyDebug(in[l] >= A->cmap->n)) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
5415: col = in[l];
5416: if (roworiented) value = v[l + k*n];
5417: else value = v[k + l*m];
5419: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
5421: if (col <= lastcol) low = 0;
5422: else high = nrow;
5423: lastcol = col;
5424: while (high-low > 5) {
5425: t = (low+high)/2;
5426: if (rp[t] > col) high = t;
5427: else low = t;
5428: }
5429: for (i=low; i<high; i++) {
5430: if (rp[i] > col) break;
5431: if (rp[i] == col) {
5432: if (is == ADD_VALUES) ap[i] += value;
5433: else ap[i] = value;
5434: goto noinsert;
5435: }
5436: }
5437: if (value == 0.0 && ignorezeroentries) goto noinsert;
5438: if (nonew == 1) goto noinsert;
5439: if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
5440: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
5441: N = nrow++ - 1; a->nz++; high++;
5442: /* shift up all the later entries in this row */
5443: for (ii=N; ii>=i; ii--) {
5444: rp[ii+1] = rp[ii];
5445: ap[ii+1] = ap[ii];
5446: }
5447: rp[i] = col;
5448: ap[i] = value;
5449: A->nonzerostate++;
5450: noinsert:;
5451: low = i + 1;
5452: }
5453: ailen[row] = nrow;
5454: }
5455: PetscFunctionReturnVoid();
5456: }