Actual source code: aij.c
petsc-3.6.4 2016-04-12
2: /*
3: Defines the basic matrix operations for the AIJ (compressed row)
4: matrix storage format.
5: */
8: #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
9: #include <petscblaslapack.h>
10: #include <petscbt.h>
11: #include <petsc/private/kernels/blocktranspose.h>
15: PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms)
16: {
18: PetscInt i,m,n;
19: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
22: MatGetSize(A,&m,&n);
23: PetscMemzero(norms,n*sizeof(PetscReal));
24: if (type == NORM_2) {
25: for (i=0; i<aij->i[m]; i++) {
26: norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]);
27: }
28: } else if (type == NORM_1) {
29: for (i=0; i<aij->i[m]; i++) {
30: norms[aij->j[i]] += PetscAbsScalar(aij->a[i]);
31: }
32: } else if (type == NORM_INFINITY) {
33: for (i=0; i<aij->i[m]; i++) {
34: norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]);
35: }
36: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType");
38: if (type == NORM_2) {
39: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
40: }
41: return(0);
42: }
46: PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is)
47: {
48: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
49: PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs;
50: const PetscInt *jj = a->j,*ii = a->i;
51: PetscInt *rows;
52: PetscErrorCode ierr;
55: for (i=0; i<m; i++) {
56: if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) {
57: cnt++;
58: }
59: }
60: PetscMalloc1(cnt,&rows);
61: cnt = 0;
62: for (i=0; i<m; i++) {
63: if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) {
64: rows[cnt] = i;
65: cnt++;
66: }
67: }
68: ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);
69: return(0);
70: }
74: PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows)
75: {
76: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
77: const MatScalar *aa = a->a;
78: PetscInt i,m=A->rmap->n,cnt = 0;
79: const PetscInt *jj = a->j,*diag;
80: PetscInt *rows;
81: PetscErrorCode ierr;
84: MatMarkDiagonal_SeqAIJ(A);
85: diag = a->diag;
86: for (i=0; i<m; i++) {
87: if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
88: cnt++;
89: }
90: }
91: PetscMalloc1(cnt,&rows);
92: cnt = 0;
93: for (i=0; i<m; i++) {
94: if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
95: rows[cnt++] = i;
96: }
97: }
98: *nrows = cnt;
99: *zrows = rows;
100: return(0);
101: }
105: PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows)
106: {
107: PetscInt nrows,*rows;
111: *zrows = NULL;
112: MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);
113: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);
114: return(0);
115: }
119: PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows)
120: {
121: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
122: const MatScalar *aa;
123: PetscInt m=A->rmap->n,cnt = 0;
124: const PetscInt *ii;
125: PetscInt n,i,j,*rows;
126: PetscErrorCode ierr;
129: *keptrows = 0;
130: ii = a->i;
131: for (i=0; i<m; i++) {
132: n = ii[i+1] - ii[i];
133: if (!n) {
134: cnt++;
135: goto ok1;
136: }
137: aa = a->a + ii[i];
138: for (j=0; j<n; j++) {
139: if (aa[j] != 0.0) goto ok1;
140: }
141: cnt++;
142: ok1:;
143: }
144: if (!cnt) return(0);
145: PetscMalloc1(A->rmap->n-cnt,&rows);
146: cnt = 0;
147: for (i=0; i<m; i++) {
148: n = ii[i+1] - ii[i];
149: if (!n) continue;
150: aa = a->a + ii[i];
151: for (j=0; j<n; j++) {
152: if (aa[j] != 0.0) {
153: rows[cnt++] = i;
154: break;
155: }
156: }
157: }
158: ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);
159: return(0);
160: }
164: PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is)
165: {
166: PetscErrorCode ierr;
167: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data;
168: PetscInt i,m = Y->rmap->n;
169: const PetscInt *diag;
170: MatScalar *aa = aij->a;
171: const PetscScalar *v;
172: PetscBool missing;
175: if (Y->assembled) {
176: MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);
177: if (!missing) {
178: diag = aij->diag;
179: VecGetArrayRead(D,&v);
180: if (is == INSERT_VALUES) {
181: for (i=0; i<m; i++) {
182: aa[diag[i]] = v[i];
183: }
184: } else {
185: for (i=0; i<m; i++) {
186: aa[diag[i]] += v[i];
187: }
188: }
189: VecRestoreArrayRead(D,&v);
190: return(0);
191: }
192: MatSeqAIJInvalidateDiagonal(Y);
193: }
194: MatDiagonalSet_Default(Y,D,is);
195: return(0);
196: }
200: PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
201: {
202: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
204: PetscInt i,ishift;
207: *m = A->rmap->n;
208: if (!ia) return(0);
209: ishift = 0;
210: if (symmetric && !A->structurally_symmetric) {
211: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);
212: } else if (oshift == 1) {
213: PetscInt *tia;
214: PetscInt nz = a->i[A->rmap->n];
215: /* malloc space and add 1 to i and j indices */
216: PetscMalloc1(A->rmap->n+1,&tia);
217: for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1;
218: *ia = tia;
219: if (ja) {
220: PetscInt *tja;
221: PetscMalloc1(nz+1,&tja);
222: for (i=0; i<nz; i++) tja[i] = a->j[i] + 1;
223: *ja = tja;
224: }
225: } else {
226: *ia = a->i;
227: if (ja) *ja = a->j;
228: }
229: return(0);
230: }
234: PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
235: {
239: if (!ia) return(0);
240: if ((symmetric && !A->structurally_symmetric) || oshift == 1) {
241: PetscFree(*ia);
242: if (ja) {PetscFree(*ja);}
243: }
244: return(0);
245: }
249: PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
250: {
251: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
253: PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n;
254: PetscInt nz = a->i[m],row,*jj,mr,col;
257: *nn = n;
258: if (!ia) return(0);
259: if (symmetric) {
260: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);
261: } else {
262: PetscCalloc1(n+1,&collengths);
263: PetscMalloc1(n+1,&cia);
264: PetscMalloc1(nz+1,&cja);
265: jj = a->j;
266: for (i=0; i<nz; i++) {
267: collengths[jj[i]]++;
268: }
269: cia[0] = oshift;
270: for (i=0; i<n; i++) {
271: cia[i+1] = cia[i] + collengths[i];
272: }
273: PetscMemzero(collengths,n*sizeof(PetscInt));
274: jj = a->j;
275: for (row=0; row<m; row++) {
276: mr = a->i[row+1] - a->i[row];
277: for (i=0; i<mr; i++) {
278: col = *jj++;
280: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
281: }
282: }
283: PetscFree(collengths);
284: *ia = cia; *ja = cja;
285: }
286: return(0);
287: }
291: PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
292: {
296: if (!ia) return(0);
298: PetscFree(*ia);
299: PetscFree(*ja);
300: return(0);
301: }
303: /*
304: MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from
305: MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output
306: spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ()
307: */
310: PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
311: {
312: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
314: PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n;
315: PetscInt nz = a->i[m],row,*jj,mr,col;
316: PetscInt *cspidx;
319: *nn = n;
320: if (!ia) return(0);
322: PetscCalloc1(n+1,&collengths);
323: PetscMalloc1(n+1,&cia);
324: PetscMalloc1(nz+1,&cja);
325: PetscMalloc1(nz+1,&cspidx);
326: jj = a->j;
327: for (i=0; i<nz; i++) {
328: collengths[jj[i]]++;
329: }
330: cia[0] = oshift;
331: for (i=0; i<n; i++) {
332: cia[i+1] = cia[i] + collengths[i];
333: }
334: PetscMemzero(collengths,n*sizeof(PetscInt));
335: jj = a->j;
336: for (row=0; row<m; row++) {
337: mr = a->i[row+1] - a->i[row];
338: for (i=0; i<mr; i++) {
339: col = *jj++;
340: cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */
341: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
342: }
343: }
344: PetscFree(collengths);
345: *ia = cia; *ja = cja;
346: *spidx = cspidx;
347: return(0);
348: }
352: PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
353: {
357: MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);
358: PetscFree(*spidx);
359: return(0);
360: }
364: PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[])
365: {
366: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
367: PetscInt *ai = a->i;
371: PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));
372: return(0);
373: }
375: /*
376: MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions
378: - a single row of values is set with each call
379: - no row or column indices are negative or (in error) larger than the number of rows or columns
380: - the values are always added to the matrix, not set
381: - no new locations are introduced in the nonzero structure of the matrix
383: This does NOT assume the global column indices are sorted
385: */
387: #include <petsc/private/isimpl.h>
390: PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
391: {
392: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
393: PetscInt low,high,t,row,nrow,i,col,l;
394: const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j;
395: PetscInt lastcol = -1;
396: MatScalar *ap,value,*aa = a->a;
397: const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices;
399: row = ridx[im[0]];
400: rp = aj + ai[row];
401: ap = aa + ai[row];
402: nrow = ailen[row];
403: low = 0;
404: high = nrow;
405: for (l=0; l<n; l++) { /* loop over added columns */
406: col = cidx[in[l]];
407: value = v[l];
409: if (col <= lastcol) low = 0;
410: else high = nrow;
411: lastcol = col;
412: while (high-low > 5) {
413: t = (low+high)/2;
414: if (rp[t] > col) high = t;
415: else low = t;
416: }
417: for (i=low; i<high; i++) {
418: if (rp[i] == col) {
419: ap[i] += value;
420: low = i + 1;
421: break;
422: }
423: }
424: }
425: return 0;
426: }
430: PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
431: {
432: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
433: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
434: PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen;
436: PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1;
437: MatScalar *ap,value,*aa = a->a;
438: PetscBool ignorezeroentries = a->ignorezeroentries;
439: PetscBool roworiented = a->roworiented;
442: for (k=0; k<m; k++) { /* loop over added rows */
443: row = im[k];
444: if (row < 0) continue;
445: #if defined(PETSC_USE_DEBUG)
446: 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);
447: #endif
448: rp = aj + ai[row]; ap = aa + ai[row];
449: rmax = imax[row]; nrow = ailen[row];
450: low = 0;
451: high = nrow;
452: for (l=0; l<n; l++) { /* loop over added columns */
453: if (in[l] < 0) continue;
454: #if defined(PETSC_USE_DEBUG)
455: 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);
456: #endif
457: col = in[l];
458: if (roworiented) {
459: value = v[l + k*n];
460: } else {
461: value = v[k + l*m];
462: }
463: if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES)) continue;
465: if (col <= lastcol) low = 0;
466: else high = nrow;
467: lastcol = col;
468: while (high-low > 5) {
469: t = (low+high)/2;
470: if (rp[t] > col) high = t;
471: else low = t;
472: }
473: for (i=low; i<high; i++) {
474: if (rp[i] > col) break;
475: if (rp[i] == col) {
476: if (is == ADD_VALUES) ap[i] += value;
477: else ap[i] = value;
478: low = i + 1;
479: goto noinsert;
480: }
481: }
482: if (value == 0.0 && ignorezeroentries) goto noinsert;
483: if (nonew == 1) goto noinsert;
484: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col);
485: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
486: N = nrow++ - 1; a->nz++; high++;
487: /* shift up all the later entries in this row */
488: for (ii=N; ii>=i; ii--) {
489: rp[ii+1] = rp[ii];
490: ap[ii+1] = ap[ii];
491: }
492: rp[i] = col;
493: ap[i] = value;
494: low = i + 1;
495: A->nonzerostate++;
496: noinsert:;
497: }
498: ailen[row] = nrow;
499: }
500: return(0);
501: }
506: PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
507: {
508: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
509: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
510: PetscInt *ai = a->i,*ailen = a->ilen;
511: MatScalar *ap,*aa = a->a;
514: for (k=0; k<m; k++) { /* loop over rows */
515: row = im[k];
516: if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
517: 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);
518: rp = aj + ai[row]; ap = aa + ai[row];
519: nrow = ailen[row];
520: for (l=0; l<n; l++) { /* loop over columns */
521: if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
522: 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);
523: col = in[l];
524: high = nrow; low = 0; /* assume unsorted */
525: while (high-low > 5) {
526: t = (low+high)/2;
527: if (rp[t] > col) high = t;
528: else low = t;
529: }
530: for (i=low; i<high; i++) {
531: if (rp[i] > col) break;
532: if (rp[i] == col) {
533: *v++ = ap[i];
534: goto finished;
535: }
536: }
537: *v++ = 0.0;
538: finished:;
539: }
540: }
541: return(0);
542: }
547: PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer)
548: {
549: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
551: PetscInt i,*col_lens;
552: int fd;
553: FILE *file;
556: PetscViewerBinaryGetDescriptor(viewer,&fd);
557: PetscMalloc1(4+A->rmap->n,&col_lens);
559: col_lens[0] = MAT_FILE_CLASSID;
560: col_lens[1] = A->rmap->n;
561: col_lens[2] = A->cmap->n;
562: col_lens[3] = a->nz;
564: /* store lengths of each row and write (including header) to file */
565: for (i=0; i<A->rmap->n; i++) {
566: col_lens[4+i] = a->i[i+1] - a->i[i];
567: }
568: PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);
569: PetscFree(col_lens);
571: /* store column indices (zero start index) */
572: PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);
574: /* store nonzero values */
575: PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);
577: PetscViewerBinaryGetInfoPointer(viewer,&file);
578: if (file) {
579: fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs));
580: }
581: return(0);
582: }
584: extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer);
588: PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer)
589: {
590: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
591: PetscErrorCode ierr;
592: PetscInt i,j,m = A->rmap->n;
593: const char *name;
594: PetscViewerFormat format;
597: PetscViewerGetFormat(viewer,&format);
598: if (format == PETSC_VIEWER_ASCII_MATLAB) {
599: PetscInt nofinalvalue = 0;
600: if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) {
601: /* Need a dummy value to ensure the dimension of the matrix. */
602: nofinalvalue = 1;
603: }
604: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
605: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);
606: PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);
607: #if defined(PETSC_USE_COMPLEX)
608: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,4);\n",a->nz+nofinalvalue);
609: #else
610: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);
611: #endif
612: PetscViewerASCIIPrintf(viewer,"zzz = [\n");
614: for (i=0; i<m; i++) {
615: for (j=a->i[i]; j<a->i[i+1]; j++) {
616: #if defined(PETSC_USE_COMPLEX)
617: 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]));
618: #else
619: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);
620: #endif
621: }
622: }
623: if (nofinalvalue) {
624: #if defined(PETSC_USE_COMPLEX)
625: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",m,A->cmap->n,0.,0.);
626: #else
627: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);
628: #endif
629: }
630: PetscObjectGetName((PetscObject)A,&name);
631: PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);
632: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
633: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
634: return(0);
635: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
636: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
637: for (i=0; i<m; i++) {
638: PetscViewerASCIIPrintf(viewer,"row %D:",i);
639: for (j=a->i[i]; j<a->i[i+1]; j++) {
640: #if defined(PETSC_USE_COMPLEX)
641: if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) {
642: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
643: } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) {
644: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));
645: } else if (PetscRealPart(a->a[j]) != 0.0) {
646: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
647: }
648: #else
649: if (a->a[j] != 0.0) {PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);}
650: #endif
651: }
652: PetscViewerASCIIPrintf(viewer,"\n");
653: }
654: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
655: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
656: PetscInt nzd=0,fshift=1,*sptr;
657: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
658: PetscMalloc1(m+1,&sptr);
659: for (i=0; i<m; i++) {
660: sptr[i] = nzd+1;
661: for (j=a->i[i]; j<a->i[i+1]; j++) {
662: if (a->j[j] >= i) {
663: #if defined(PETSC_USE_COMPLEX)
664: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++;
665: #else
666: if (a->a[j] != 0.0) nzd++;
667: #endif
668: }
669: }
670: }
671: sptr[m] = nzd+1;
672: PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);
673: for (i=0; i<m+1; i+=6) {
674: if (i+4<m) {
675: 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]);
676: } else if (i+3<m) {
677: PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);
678: } else if (i+2<m) {
679: PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);
680: } else if (i+1<m) {
681: PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);
682: } else if (i<m) {
683: PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);
684: } else {
685: PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);
686: }
687: }
688: PetscViewerASCIIPrintf(viewer,"\n");
689: PetscFree(sptr);
690: for (i=0; i<m; i++) {
691: for (j=a->i[i]; j<a->i[i+1]; j++) {
692: if (a->j[j] >= i) {PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);}
693: }
694: PetscViewerASCIIPrintf(viewer,"\n");
695: }
696: PetscViewerASCIIPrintf(viewer,"\n");
697: for (i=0; i<m; i++) {
698: for (j=a->i[i]; j<a->i[i+1]; j++) {
699: if (a->j[j] >= i) {
700: #if defined(PETSC_USE_COMPLEX)
701: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) {
702: PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
703: }
704: #else
705: if (a->a[j] != 0.0) {PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);}
706: #endif
707: }
708: }
709: PetscViewerASCIIPrintf(viewer,"\n");
710: }
711: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
712: } else if (format == PETSC_VIEWER_ASCII_DENSE) {
713: PetscInt cnt = 0,jcnt;
714: PetscScalar value;
715: #if defined(PETSC_USE_COMPLEX)
716: PetscBool realonly = PETSC_TRUE;
718: for (i=0; i<a->i[m]; i++) {
719: if (PetscImaginaryPart(a->a[i]) != 0.0) {
720: realonly = PETSC_FALSE;
721: break;
722: }
723: }
724: #endif
726: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
727: for (i=0; i<m; i++) {
728: jcnt = 0;
729: for (j=0; j<A->cmap->n; j++) {
730: if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) {
731: value = a->a[cnt++];
732: jcnt++;
733: } else {
734: value = 0.0;
735: }
736: #if defined(PETSC_USE_COMPLEX)
737: if (realonly) {
738: PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));
739: } else {
740: PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));
741: }
742: #else
743: PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);
744: #endif
745: }
746: PetscViewerASCIIPrintf(viewer,"\n");
747: }
748: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
749: } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) {
750: PetscInt fshift=1;
751: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
752: #if defined(PETSC_USE_COMPLEX)
753: PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n");
754: #else
755: PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n");
756: #endif
757: PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);
758: for (i=0; i<m; i++) {
759: for (j=a->i[i]; j<a->i[i+1]; j++) {
760: #if defined(PETSC_USE_COMPLEX)
761: PetscViewerASCIIPrintf(viewer,"%D %D %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
762: #else
763: PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);
764: #endif
765: }
766: }
767: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
768: } else {
769: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
770: if (A->factortype) {
771: for (i=0; i<m; i++) {
772: PetscViewerASCIIPrintf(viewer,"row %D:",i);
773: /* L part */
774: for (j=a->i[i]; j<a->i[i+1]; j++) {
775: #if defined(PETSC_USE_COMPLEX)
776: if (PetscImaginaryPart(a->a[j]) > 0.0) {
777: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
778: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
779: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));
780: } else {
781: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
782: }
783: #else
784: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);
785: #endif
786: }
787: /* diagonal */
788: j = a->diag[i];
789: #if defined(PETSC_USE_COMPLEX)
790: if (PetscImaginaryPart(a->a[j]) > 0.0) {
791: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]));
792: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
793: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])));
794: } else {
795: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));
796: }
797: #else
798: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)(1.0/a->a[j]));
799: #endif
801: /* U part */
802: for (j=a->diag[i+1]+1; j<a->diag[i]; j++) {
803: #if defined(PETSC_USE_COMPLEX)
804: if (PetscImaginaryPart(a->a[j]) > 0.0) {
805: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
806: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
807: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));
808: } else {
809: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
810: }
811: #else
812: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);
813: #endif
814: }
815: PetscViewerASCIIPrintf(viewer,"\n");
816: }
817: } else {
818: for (i=0; i<m; i++) {
819: PetscViewerASCIIPrintf(viewer,"row %D:",i);
820: for (j=a->i[i]; j<a->i[i+1]; j++) {
821: #if defined(PETSC_USE_COMPLEX)
822: if (PetscImaginaryPart(a->a[j]) > 0.0) {
823: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));
824: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
825: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));
826: } else {
827: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));
828: }
829: #else
830: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);
831: #endif
832: }
833: PetscViewerASCIIPrintf(viewer,"\n");
834: }
835: }
836: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
837: }
838: PetscViewerFlush(viewer);
839: return(0);
840: }
842: #include <petscdraw.h>
845: PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
846: {
847: Mat A = (Mat) Aa;
848: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
849: PetscErrorCode ierr;
850: PetscInt i,j,m = A->rmap->n,color;
851: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0;
852: PetscViewer viewer;
853: PetscViewerFormat format;
856: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
857: PetscViewerGetFormat(viewer,&format);
859: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
860: /* loop over matrix elements drawing boxes */
862: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
863: /* Blue for negative, Cyan for zero and Red for positive */
864: color = PETSC_DRAW_BLUE;
865: for (i=0; i<m; i++) {
866: y_l = m - i - 1.0; y_r = y_l + 1.0;
867: for (j=a->i[i]; j<a->i[i+1]; j++) {
868: x_l = a->j[j]; x_r = x_l + 1.0;
869: if (PetscRealPart(a->a[j]) >= 0.) continue;
870: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
871: }
872: }
873: color = PETSC_DRAW_CYAN;
874: for (i=0; i<m; i++) {
875: y_l = m - i - 1.0; y_r = y_l + 1.0;
876: for (j=a->i[i]; j<a->i[i+1]; j++) {
877: x_l = a->j[j]; x_r = x_l + 1.0;
878: if (a->a[j] != 0.) continue;
879: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
880: }
881: }
882: color = PETSC_DRAW_RED;
883: for (i=0; i<m; i++) {
884: y_l = m - i - 1.0; y_r = y_l + 1.0;
885: for (j=a->i[i]; j<a->i[i+1]; j++) {
886: x_l = a->j[j]; x_r = x_l + 1.0;
887: if (PetscRealPart(a->a[j]) <= 0.) continue;
888: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
889: }
890: }
891: } else {
892: /* use contour shading to indicate magnitude of values */
893: /* first determine max of all nonzero values */
894: PetscInt nz = a->nz,count;
895: PetscDraw popup;
896: PetscReal scale;
898: for (i=0; i<nz; i++) {
899: if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
900: }
901: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
902: PetscDrawGetPopup(draw,&popup);
903: if (popup) {
904: PetscDrawScalePopup(popup,0.0,maxv);
905: }
906: count = 0;
907: for (i=0; i<m; i++) {
908: y_l = m - i - 1.0; y_r = y_l + 1.0;
909: for (j=a->i[i]; j<a->i[i+1]; j++) {
910: x_l = a->j[j]; x_r = x_l + 1.0;
911: color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count]));
912: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
913: count++;
914: }
915: }
916: }
917: return(0);
918: }
920: #include <petscdraw.h>
923: PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer)
924: {
926: PetscDraw draw;
927: PetscReal xr,yr,xl,yl,h,w;
928: PetscBool isnull;
931: PetscViewerDrawGetDraw(viewer,0,&draw);
932: PetscDrawIsNull(draw,&isnull);
933: if (isnull) return(0);
935: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
936: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
937: xr += w; yr += h; xl = -w; yl = -h;
938: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
939: PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);
940: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
941: return(0);
942: }
946: PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer)
947: {
949: PetscBool iascii,isbinary,isdraw;
952: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
953: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
954: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
955: if (iascii) {
956: MatView_SeqAIJ_ASCII(A,viewer);
957: } else if (isbinary) {
958: MatView_SeqAIJ_Binary(A,viewer);
959: } else if (isdraw) {
960: MatView_SeqAIJ_Draw(A,viewer);
961: }
962: MatView_SeqAIJ_Inode(A,viewer);
963: return(0);
964: }
968: PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode)
969: {
970: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
972: PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
973: PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0;
974: MatScalar *aa = a->a,*ap;
975: PetscReal ratio = 0.6;
978: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
980: if (m) rmax = ailen[0]; /* determine row with most nonzeros */
981: for (i=1; i<m; i++) {
982: /* move each row back by the amount of empty slots (fshift) before it*/
983: fshift += imax[i-1] - ailen[i-1];
984: rmax = PetscMax(rmax,ailen[i]);
985: if (fshift) {
986: ip = aj + ai[i];
987: ap = aa + ai[i];
988: N = ailen[i];
989: for (j=0; j<N; j++) {
990: ip[j-fshift] = ip[j];
991: ap[j-fshift] = ap[j];
992: }
993: }
994: ai[i] = ai[i-1] + ailen[i-1];
995: }
996: if (m) {
997: fshift += imax[m-1] - ailen[m-1];
998: ai[m] = ai[m-1] + ailen[m-1];
999: }
1001: /* reset ilen and imax for each row */
1002: a->nonzerorowcnt = 0;
1003: for (i=0; i<m; i++) {
1004: ailen[i] = imax[i] = ai[i+1] - ai[i];
1005: a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0);
1006: }
1007: a->nz = ai[m];
1008: 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);
1010: MatMarkDiagonal_SeqAIJ(A);
1011: PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);
1012: PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);
1013: PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);
1015: A->info.mallocs += a->reallocs;
1016: a->reallocs = 0;
1017: A->info.nz_unneeded = (PetscReal)fshift;
1018: a->rmax = rmax;
1020: MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);
1021: MatAssemblyEnd_SeqAIJ_Inode(A,mode);
1022: MatSeqAIJInvalidateDiagonal(A);
1023: return(0);
1024: }
1028: PetscErrorCode MatRealPart_SeqAIJ(Mat A)
1029: {
1030: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1031: PetscInt i,nz = a->nz;
1032: MatScalar *aa = a->a;
1036: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1037: MatSeqAIJInvalidateDiagonal(A);
1038: return(0);
1039: }
1043: PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A)
1044: {
1045: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1046: PetscInt i,nz = a->nz;
1047: MatScalar *aa = a->a;
1051: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1052: MatSeqAIJInvalidateDiagonal(A);
1053: return(0);
1054: }
1058: PetscErrorCode MatZeroEntries_SeqAIJ(Mat A)
1059: {
1060: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1064: PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));
1065: MatSeqAIJInvalidateDiagonal(A);
1066: return(0);
1067: }
1071: PetscErrorCode MatDestroy_SeqAIJ(Mat A)
1072: {
1073: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1077: #if defined(PETSC_USE_LOG)
1078: PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz);
1079: #endif
1080: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
1081: ISDestroy(&a->row);
1082: ISDestroy(&a->col);
1083: PetscFree(a->diag);
1084: PetscFree(a->ibdiag);
1085: PetscFree2(a->imax,a->ilen);
1086: PetscFree3(a->idiag,a->mdiag,a->ssor_work);
1087: PetscFree(a->solve_work);
1088: ISDestroy(&a->icol);
1089: PetscFree(a->saved_values);
1090: ISColoringDestroy(&a->coloring);
1091: PetscFree2(a->compressedrow.i,a->compressedrow.rindex);
1092: PetscFree(a->matmult_abdense);
1094: MatDestroy_SeqAIJ_Inode(A);
1095: PetscFree(A->data);
1097: PetscObjectChangeTypeName((PetscObject)A,0);
1098: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);
1099: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
1100: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
1101: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);
1102: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);
1103: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);
1104: #if defined(PETSC_HAVE_ELEMENTAL)
1105: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL);
1106: #endif
1107: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL);
1108: PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);
1109: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);
1110: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);
1111: PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);
1112: return(0);
1113: }
1117: PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg)
1118: {
1119: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1123: switch (op) {
1124: case MAT_ROW_ORIENTED:
1125: a->roworiented = flg;
1126: break;
1127: case MAT_KEEP_NONZERO_PATTERN:
1128: a->keepnonzeropattern = flg;
1129: break;
1130: case MAT_NEW_NONZERO_LOCATIONS:
1131: a->nonew = (flg ? 0 : 1);
1132: break;
1133: case MAT_NEW_NONZERO_LOCATION_ERR:
1134: a->nonew = (flg ? -1 : 0);
1135: break;
1136: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1137: a->nonew = (flg ? -2 : 0);
1138: break;
1139: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1140: a->nounused = (flg ? -1 : 0);
1141: break;
1142: case MAT_IGNORE_ZERO_ENTRIES:
1143: a->ignorezeroentries = flg;
1144: break;
1145: case MAT_SPD:
1146: case MAT_SYMMETRIC:
1147: case MAT_STRUCTURALLY_SYMMETRIC:
1148: case MAT_HERMITIAN:
1149: case MAT_SYMMETRY_ETERNAL:
1150: /* These options are handled directly by MatSetOption() */
1151: break;
1152: case MAT_NEW_DIAGONALS:
1153: case MAT_IGNORE_OFF_PROC_ENTRIES:
1154: case MAT_USE_HASH_TABLE:
1155: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1156: break;
1157: case MAT_USE_INODES:
1158: /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */
1159: break;
1160: default:
1161: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
1162: }
1163: MatSetOption_SeqAIJ_Inode(A,op,flg);
1164: return(0);
1165: }
1169: PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v)
1170: {
1171: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1173: PetscInt i,j,n,*ai=a->i,*aj=a->j,nz;
1174: PetscScalar *aa=a->a,*x,zero=0.0;
1177: VecGetLocalSize(v,&n);
1178: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1180: if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) {
1181: PetscInt *diag=a->diag;
1182: VecGetArray(v,&x);
1183: for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]];
1184: VecRestoreArray(v,&x);
1185: return(0);
1186: }
1188: VecSet(v,zero);
1189: VecGetArray(v,&x);
1190: for (i=0; i<n; i++) {
1191: nz = ai[i+1] - ai[i];
1192: if (!nz) x[i] = 0.0;
1193: for (j=ai[i]; j<ai[i+1]; j++) {
1194: if (aj[j] == i) {
1195: x[i] = aa[j];
1196: break;
1197: }
1198: }
1199: }
1200: VecRestoreArray(v,&x);
1201: return(0);
1202: }
1204: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1207: PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy)
1208: {
1209: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1210: PetscScalar *y;
1211: const PetscScalar *x;
1212: PetscErrorCode ierr;
1213: PetscInt m = A->rmap->n;
1214: #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1215: const MatScalar *v;
1216: PetscScalar alpha;
1217: PetscInt n,i,j;
1218: const PetscInt *idx,*ii,*ridx=NULL;
1219: Mat_CompressedRow cprow = a->compressedrow;
1220: PetscBool usecprow = cprow.use;
1221: #endif
1224: if (zz != yy) {VecCopy(zz,yy);}
1225: VecGetArrayRead(xx,&x);
1226: VecGetArray(yy,&y);
1228: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1229: fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y);
1230: #else
1231: if (usecprow) {
1232: m = cprow.nrows;
1233: ii = cprow.i;
1234: ridx = cprow.rindex;
1235: } else {
1236: ii = a->i;
1237: }
1238: for (i=0; i<m; i++) {
1239: idx = a->j + ii[i];
1240: v = a->a + ii[i];
1241: n = ii[i+1] - ii[i];
1242: if (usecprow) {
1243: alpha = x[ridx[i]];
1244: } else {
1245: alpha = x[i];
1246: }
1247: for (j=0; j<n; j++) y[idx[j]] += alpha*v[j];
1248: }
1249: #endif
1250: PetscLogFlops(2.0*a->nz);
1251: VecRestoreArrayRead(xx,&x);
1252: VecRestoreArray(yy,&y);
1253: return(0);
1254: }
1258: PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy)
1259: {
1263: VecSet(yy,0.0);
1264: MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);
1265: return(0);
1266: }
1268: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1272: PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy)
1273: {
1274: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1275: PetscScalar *y;
1276: const PetscScalar *x;
1277: const MatScalar *aa;
1278: PetscErrorCode ierr;
1279: PetscInt m=A->rmap->n;
1280: const PetscInt *aj,*ii,*ridx=NULL;
1281: PetscInt n,i;
1282: PetscScalar sum;
1283: PetscBool usecprow=a->compressedrow.use;
1285: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1286: #pragma disjoint(*x,*y,*aa)
1287: #endif
1290: VecGetArrayRead(xx,&x);
1291: VecGetArray(yy,&y);
1292: aj = a->j;
1293: aa = a->a;
1294: ii = a->i;
1295: if (usecprow) { /* use compressed row format */
1296: PetscMemzero(y,m*sizeof(PetscScalar));
1297: m = a->compressedrow.nrows;
1298: ii = a->compressedrow.i;
1299: ridx = a->compressedrow.rindex;
1300: for (i=0; i<m; i++) {
1301: n = ii[i+1] - ii[i];
1302: aj = a->j + ii[i];
1303: aa = a->a + ii[i];
1304: sum = 0.0;
1305: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1306: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1307: y[*ridx++] = sum;
1308: }
1309: } else { /* do not use compressed row format */
1310: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ)
1311: fortranmultaij_(&m,x,ii,aj,aa,y);
1312: #else
1313: for (i=0; i<m; i++) {
1314: n = ii[i+1] - ii[i];
1315: aj = a->j + ii[i];
1316: aa = a->a + ii[i];
1317: sum = 0.0;
1318: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1319: y[i] = sum;
1320: }
1321: #endif
1322: }
1323: PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);
1324: VecRestoreArrayRead(xx,&x);
1325: VecRestoreArray(yy,&y);
1326: return(0);
1327: }
1331: PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy)
1332: {
1333: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1334: PetscScalar *y;
1335: const PetscScalar *x;
1336: const MatScalar *aa;
1337: PetscErrorCode ierr;
1338: PetscInt m=A->rmap->n;
1339: const PetscInt *aj,*ii,*ridx=NULL;
1340: PetscInt n,i,nonzerorow=0;
1341: PetscScalar sum;
1342: PetscBool usecprow=a->compressedrow.use;
1344: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1345: #pragma disjoint(*x,*y,*aa)
1346: #endif
1349: VecGetArrayRead(xx,&x);
1350: VecGetArray(yy,&y);
1351: aj = a->j;
1352: aa = a->a;
1353: ii = a->i;
1354: if (usecprow) { /* use compressed row format */
1355: m = a->compressedrow.nrows;
1356: ii = a->compressedrow.i;
1357: ridx = a->compressedrow.rindex;
1358: for (i=0; i<m; i++) {
1359: n = ii[i+1] - ii[i];
1360: aj = a->j + ii[i];
1361: aa = a->a + ii[i];
1362: sum = 0.0;
1363: nonzerorow += (n>0);
1364: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1365: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1366: y[*ridx++] = sum;
1367: }
1368: } else { /* do not use compressed row format */
1369: for (i=0; i<m; i++) {
1370: n = ii[i+1] - ii[i];
1371: aj = a->j + ii[i];
1372: aa = a->a + ii[i];
1373: sum = 0.0;
1374: nonzerorow += (n>0);
1375: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1376: y[i] = sum;
1377: }
1378: }
1379: PetscLogFlops(2.0*a->nz - nonzerorow);
1380: VecRestoreArrayRead(xx,&x);
1381: VecRestoreArray(yy,&y);
1382: return(0);
1383: }
1387: PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz)
1388: {
1389: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1390: PetscScalar *y,*z;
1391: const PetscScalar *x;
1392: const MatScalar *aa;
1393: PetscErrorCode ierr;
1394: PetscInt m = A->rmap->n,*aj,*ii;
1395: PetscInt n,i,*ridx=NULL;
1396: PetscScalar sum;
1397: PetscBool usecprow=a->compressedrow.use;
1400: VecGetArrayRead(xx,&x);
1401: VecGetArrayPair(yy,zz,&y,&z);
1403: aj = a->j;
1404: aa = a->a;
1405: ii = a->i;
1406: if (usecprow) { /* use compressed row format */
1407: if (zz != yy) {
1408: PetscMemcpy(z,y,m*sizeof(PetscScalar));
1409: }
1410: m = a->compressedrow.nrows;
1411: ii = a->compressedrow.i;
1412: ridx = a->compressedrow.rindex;
1413: for (i=0; i<m; i++) {
1414: n = ii[i+1] - ii[i];
1415: aj = a->j + ii[i];
1416: aa = a->a + ii[i];
1417: sum = y[*ridx];
1418: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1419: z[*ridx++] = sum;
1420: }
1421: } else { /* do not use compressed row format */
1422: for (i=0; i<m; i++) {
1423: n = ii[i+1] - ii[i];
1424: aj = a->j + ii[i];
1425: aa = a->a + ii[i];
1426: sum = y[i];
1427: PetscSparseDenseMaxDot(sum,x,aa,aj,n);
1428: z[i] = sum;
1429: }
1430: }
1431: PetscLogFlops(2.0*a->nz);
1432: VecRestoreArrayRead(xx,&x);
1433: VecRestoreArrayPair(yy,zz,&y,&z);
1434: return(0);
1435: }
1437: #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h>
1440: PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz)
1441: {
1442: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1443: PetscScalar *y,*z;
1444: const PetscScalar *x;
1445: const MatScalar *aa;
1446: PetscErrorCode ierr;
1447: const PetscInt *aj,*ii,*ridx=NULL;
1448: PetscInt m = A->rmap->n,n,i;
1449: PetscScalar sum;
1450: PetscBool usecprow=a->compressedrow.use;
1453: VecGetArrayRead(xx,&x);
1454: VecGetArrayPair(yy,zz,&y,&z);
1456: aj = a->j;
1457: aa = a->a;
1458: ii = a->i;
1459: if (usecprow) { /* use compressed row format */
1460: if (zz != yy) {
1461: PetscMemcpy(z,y,m*sizeof(PetscScalar));
1462: }
1463: m = a->compressedrow.nrows;
1464: ii = a->compressedrow.i;
1465: ridx = a->compressedrow.rindex;
1466: for (i=0; i<m; i++) {
1467: n = ii[i+1] - ii[i];
1468: aj = a->j + ii[i];
1469: aa = a->a + ii[i];
1470: sum = y[*ridx];
1471: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1472: z[*ridx++] = sum;
1473: }
1474: } else { /* do not use compressed row format */
1475: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ)
1476: fortranmultaddaij_(&m,x,ii,aj,aa,y,z);
1477: #else
1478: for (i=0; i<m; i++) {
1479: n = ii[i+1] - ii[i];
1480: aj = a->j + ii[i];
1481: aa = a->a + ii[i];
1482: sum = y[i];
1483: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1484: z[i] = sum;
1485: }
1486: #endif
1487: }
1488: PetscLogFlops(2.0*a->nz);
1489: VecRestoreArrayRead(xx,&x);
1490: VecRestoreArrayPair(yy,zz,&y,&z);
1491: #if defined(PETSC_HAVE_CUSP)
1492: /*
1493: VecView(xx,0);
1494: VecView(zz,0);
1495: MatView(A,0);
1496: */
1497: #endif
1498: return(0);
1499: }
1501: /*
1502: Adds diagonal pointers to sparse matrix structure.
1503: */
1506: PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A)
1507: {
1508: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1510: PetscInt i,j,m = A->rmap->n;
1513: if (!a->diag) {
1514: PetscMalloc1(m,&a->diag);
1515: PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));
1516: }
1517: for (i=0; i<A->rmap->n; i++) {
1518: a->diag[i] = a->i[i+1];
1519: for (j=a->i[i]; j<a->i[i+1]; j++) {
1520: if (a->j[j] == i) {
1521: a->diag[i] = j;
1522: break;
1523: }
1524: }
1525: }
1526: return(0);
1527: }
1529: /*
1530: Checks for missing diagonals
1531: */
1534: PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d)
1535: {
1536: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1537: PetscInt *diag,*ii = a->i,i;
1540: *missing = PETSC_FALSE;
1541: if (A->rmap->n > 0 && !ii) {
1542: *missing = PETSC_TRUE;
1543: if (d) *d = 0;
1544: PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
1545: } else {
1546: diag = a->diag;
1547: for (i=0; i<A->rmap->n; i++) {
1548: if (diag[i] >= ii[i+1]) {
1549: *missing = PETSC_TRUE;
1550: if (d) *d = i;
1551: PetscInfo1(A,"Matrix is missing diagonal number %D\n",i);
1552: break;
1553: }
1554: }
1555: }
1556: return(0);
1557: }
1561: /*
1562: Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways
1563: */
1564: PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift)
1565: {
1566: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
1568: PetscInt i,*diag,m = A->rmap->n;
1569: MatScalar *v = a->a;
1570: PetscScalar *idiag,*mdiag;
1573: if (a->idiagvalid) return(0);
1574: MatMarkDiagonal_SeqAIJ(A);
1575: diag = a->diag;
1576: if (!a->idiag) {
1577: PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);
1578: PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));
1579: v = a->a;
1580: }
1581: mdiag = a->mdiag;
1582: idiag = a->idiag;
1584: if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) {
1585: for (i=0; i<m; i++) {
1586: mdiag[i] = v[diag[i]];
1587: if (!PetscAbsScalar(mdiag[i]) && !PetscRealPart(fshift)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i);
1588: idiag[i] = 1.0/v[diag[i]];
1589: }
1590: PetscLogFlops(m);
1591: } else {
1592: for (i=0; i<m; i++) {
1593: mdiag[i] = v[diag[i]];
1594: idiag[i] = omega/(fshift + v[diag[i]]);
1595: }
1596: PetscLogFlops(2.0*m);
1597: }
1598: a->idiagvalid = PETSC_TRUE;
1599: return(0);
1600: }
1602: #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h>
1605: PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1606: {
1607: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1608: PetscScalar *x,d,sum,*t,scale;
1609: const MatScalar *v = a->a,*idiag=0,*mdiag;
1610: const PetscScalar *b, *bs,*xb, *ts;
1611: PetscErrorCode ierr;
1612: PetscInt n = A->cmap->n,m = A->rmap->n,i;
1613: const PetscInt *idx,*diag;
1616: its = its*lits;
1618: if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */
1619: if (!a->idiagvalid) {MatInvertDiagonal_SeqAIJ(A,omega,fshift);}
1620: a->fshift = fshift;
1621: a->omega = omega;
1623: diag = a->diag;
1624: t = a->ssor_work;
1625: idiag = a->idiag;
1626: mdiag = a->mdiag;
1628: VecGetArray(xx,&x);
1629: VecGetArrayRead(bb,&b);
1630: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
1631: if (flag == SOR_APPLY_UPPER) {
1632: /* apply (U + D/omega) to the vector */
1633: bs = b;
1634: for (i=0; i<m; i++) {
1635: d = fshift + mdiag[i];
1636: n = a->i[i+1] - diag[i] - 1;
1637: idx = a->j + diag[i] + 1;
1638: v = a->a + diag[i] + 1;
1639: sum = b[i]*d/omega;
1640: PetscSparseDensePlusDot(sum,bs,v,idx,n);
1641: x[i] = sum;
1642: }
1643: VecRestoreArray(xx,&x);
1644: VecRestoreArrayRead(bb,&b);
1645: PetscLogFlops(a->nz);
1646: return(0);
1647: }
1649: if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented");
1650: else if (flag & SOR_EISENSTAT) {
1651: /* Let A = L + U + D; where L is lower trianglar,
1652: U is upper triangular, E = D/omega; This routine applies
1654: (L + E)^{-1} A (U + E)^{-1}
1656: to a vector efficiently using Eisenstat's trick.
1657: */
1658: scale = (2.0/omega) - 1.0;
1660: /* x = (E + U)^{-1} b */
1661: for (i=m-1; i>=0; i--) {
1662: n = a->i[i+1] - diag[i] - 1;
1663: idx = a->j + diag[i] + 1;
1664: v = a->a + diag[i] + 1;
1665: sum = b[i];
1666: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1667: x[i] = sum*idiag[i];
1668: }
1670: /* t = b - (2*E - D)x */
1671: v = a->a;
1672: for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i];
1674: /* t = (E + L)^{-1}t */
1675: ts = t;
1676: diag = a->diag;
1677: for (i=0; i<m; i++) {
1678: n = diag[i] - a->i[i];
1679: idx = a->j + a->i[i];
1680: v = a->a + a->i[i];
1681: sum = t[i];
1682: PetscSparseDenseMinusDot(sum,ts,v,idx,n);
1683: t[i] = sum*idiag[i];
1684: /* x = x + t */
1685: x[i] += t[i];
1686: }
1688: PetscLogFlops(6.0*m-1 + 2.0*a->nz);
1689: VecRestoreArray(xx,&x);
1690: VecRestoreArrayRead(bb,&b);
1691: return(0);
1692: }
1693: if (flag & SOR_ZERO_INITIAL_GUESS) {
1694: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1695: for (i=0; i<m; i++) {
1696: n = diag[i] - a->i[i];
1697: idx = a->j + a->i[i];
1698: v = a->a + a->i[i];
1699: sum = b[i];
1700: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1701: t[i] = sum;
1702: x[i] = sum*idiag[i];
1703: }
1704: xb = t;
1705: PetscLogFlops(a->nz);
1706: } else xb = b;
1707: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1708: for (i=m-1; i>=0; i--) {
1709: n = a->i[i+1] - diag[i] - 1;
1710: idx = a->j + diag[i] + 1;
1711: v = a->a + diag[i] + 1;
1712: sum = xb[i];
1713: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1714: if (xb == b) {
1715: x[i] = sum*idiag[i];
1716: } else {
1717: x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */
1718: }
1719: }
1720: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
1721: }
1722: its--;
1723: }
1724: while (its--) {
1725: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1726: for (i=0; i<m; i++) {
1727: /* lower */
1728: n = diag[i] - a->i[i];
1729: idx = a->j + a->i[i];
1730: v = a->a + a->i[i];
1731: sum = b[i];
1732: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1733: t[i] = sum; /* save application of the lower-triangular part */
1734: /* upper */
1735: n = a->i[i+1] - diag[i] - 1;
1736: idx = a->j + diag[i] + 1;
1737: v = a->a + diag[i] + 1;
1738: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1739: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
1740: }
1741: xb = t;
1742: PetscLogFlops(2.0*a->nz);
1743: } else xb = b;
1744: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1745: for (i=m-1; i>=0; i--) {
1746: sum = xb[i];
1747: if (xb == b) {
1748: /* whole matrix (no checkpointing available) */
1749: n = a->i[i+1] - a->i[i];
1750: idx = a->j + a->i[i];
1751: v = a->a + a->i[i];
1752: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1753: x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i];
1754: } else { /* lower-triangular part has been saved, so only apply upper-triangular */
1755: n = a->i[i+1] - diag[i] - 1;
1756: idx = a->j + diag[i] + 1;
1757: v = a->a + diag[i] + 1;
1758: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1759: x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */
1760: }
1761: }
1762: if (xb == b) {
1763: PetscLogFlops(2.0*a->nz);
1764: } else {
1765: PetscLogFlops(a->nz); /* assumes 1/2 in upper */
1766: }
1767: }
1768: }
1769: VecRestoreArray(xx,&x);
1770: VecRestoreArrayRead(bb,&b);
1771: return(0);
1772: }
1777: PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info)
1778: {
1779: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1782: info->block_size = 1.0;
1783: info->nz_allocated = (double)a->maxnz;
1784: info->nz_used = (double)a->nz;
1785: info->nz_unneeded = (double)(a->maxnz - a->nz);
1786: info->assemblies = (double)A->num_ass;
1787: info->mallocs = (double)A->info.mallocs;
1788: info->memory = ((PetscObject)A)->mem;
1789: if (A->factortype) {
1790: info->fill_ratio_given = A->info.fill_ratio_given;
1791: info->fill_ratio_needed = A->info.fill_ratio_needed;
1792: info->factor_mallocs = A->info.factor_mallocs;
1793: } else {
1794: info->fill_ratio_given = 0;
1795: info->fill_ratio_needed = 0;
1796: info->factor_mallocs = 0;
1797: }
1798: return(0);
1799: }
1803: PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1804: {
1805: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1806: PetscInt i,m = A->rmap->n - 1,d = 0;
1807: PetscErrorCode ierr;
1808: const PetscScalar *xx;
1809: PetscScalar *bb;
1810: PetscBool missing;
1813: if (x && b) {
1814: VecGetArrayRead(x,&xx);
1815: VecGetArray(b,&bb);
1816: for (i=0; i<N; i++) {
1817: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1818: bb[rows[i]] = diag*xx[rows[i]];
1819: }
1820: VecRestoreArrayRead(x,&xx);
1821: VecRestoreArray(b,&bb);
1822: }
1824: if (a->keepnonzeropattern) {
1825: for (i=0; i<N; i++) {
1826: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1827: PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));
1828: }
1829: if (diag != 0.0) {
1830: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
1831: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d);
1832: for (i=0; i<N; i++) {
1833: a->a[a->diag[rows[i]]] = diag;
1834: }
1835: }
1836: } else {
1837: if (diag != 0.0) {
1838: for (i=0; i<N; i++) {
1839: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1840: if (a->ilen[rows[i]] > 0) {
1841: a->ilen[rows[i]] = 1;
1842: a->a[a->i[rows[i]]] = diag;
1843: a->j[a->i[rows[i]]] = rows[i];
1844: } else { /* in case row was completely empty */
1845: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
1846: }
1847: }
1848: } else {
1849: for (i=0; i<N; i++) {
1850: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1851: a->ilen[rows[i]] = 0;
1852: }
1853: }
1854: A->nonzerostate++;
1855: }
1856: MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);
1857: return(0);
1858: }
1862: PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1863: {
1864: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1865: PetscInt i,j,m = A->rmap->n - 1,d = 0;
1866: PetscErrorCode ierr;
1867: PetscBool missing,*zeroed,vecs = PETSC_FALSE;
1868: const PetscScalar *xx;
1869: PetscScalar *bb;
1872: if (x && b) {
1873: VecGetArrayRead(x,&xx);
1874: VecGetArray(b,&bb);
1875: vecs = PETSC_TRUE;
1876: }
1877: PetscCalloc1(A->rmap->n,&zeroed);
1878: for (i=0; i<N; i++) {
1879: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1880: PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));
1882: zeroed[rows[i]] = PETSC_TRUE;
1883: }
1884: for (i=0; i<A->rmap->n; i++) {
1885: if (!zeroed[i]) {
1886: for (j=a->i[i]; j<a->i[i+1]; j++) {
1887: if (zeroed[a->j[j]]) {
1888: if (vecs) bb[i] -= a->a[j]*xx[a->j[j]];
1889: a->a[j] = 0.0;
1890: }
1891: }
1892: } else if (vecs) bb[i] = diag*xx[i];
1893: }
1894: if (x && b) {
1895: VecRestoreArrayRead(x,&xx);
1896: VecRestoreArray(b,&bb);
1897: }
1898: PetscFree(zeroed);
1899: if (diag != 0.0) {
1900: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
1901: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d);
1902: for (i=0; i<N; i++) {
1903: a->a[a->diag[rows[i]]] = diag;
1904: }
1905: }
1906: MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);
1907: return(0);
1908: }
1912: PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1913: {
1914: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1915: PetscInt *itmp;
1918: if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row);
1920: *nz = a->i[row+1] - a->i[row];
1921: if (v) *v = a->a + a->i[row];
1922: if (idx) {
1923: itmp = a->j + a->i[row];
1924: if (*nz) *idx = itmp;
1925: else *idx = 0;
1926: }
1927: return(0);
1928: }
1930: /* remove this function? */
1933: PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1934: {
1936: return(0);
1937: }
1941: PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm)
1942: {
1943: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1944: MatScalar *v = a->a;
1945: PetscReal sum = 0.0;
1947: PetscInt i,j;
1950: if (type == NORM_FROBENIUS) {
1951: for (i=0; i<a->nz; i++) {
1952: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1953: }
1954: *nrm = PetscSqrtReal(sum);
1955: } else if (type == NORM_1) {
1956: PetscReal *tmp;
1957: PetscInt *jj = a->j;
1958: PetscCalloc1(A->cmap->n+1,&tmp);
1959: *nrm = 0.0;
1960: for (j=0; j<a->nz; j++) {
1961: tmp[*jj++] += PetscAbsScalar(*v); v++;
1962: }
1963: for (j=0; j<A->cmap->n; j++) {
1964: if (tmp[j] > *nrm) *nrm = tmp[j];
1965: }
1966: PetscFree(tmp);
1967: } else if (type == NORM_INFINITY) {
1968: *nrm = 0.0;
1969: for (j=0; j<A->rmap->n; j++) {
1970: v = a->a + a->i[j];
1971: sum = 0.0;
1972: for (i=0; i<a->i[j+1]-a->i[j]; i++) {
1973: sum += PetscAbsScalar(*v); v++;
1974: }
1975: if (sum > *nrm) *nrm = sum;
1976: }
1977: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm");
1978: return(0);
1979: }
1981: /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */
1984: PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B)
1985: {
1987: PetscInt i,j,anzj;
1988: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b;
1989: PetscInt an=A->cmap->N,am=A->rmap->N;
1990: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
1993: /* Allocate space for symbolic transpose info and work array */
1994: PetscCalloc1(an+1,&ati);
1995: PetscMalloc1(ai[am],&atj);
1996: PetscMalloc1(an,&atfill);
1998: /* Walk through aj and count ## of non-zeros in each row of A^T. */
1999: /* Note: offset by 1 for fast conversion into csr format. */
2000: for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1;
2001: /* Form ati for csr format of A^T. */
2002: for (i=0;i<an;i++) ati[i+1] += ati[i];
2004: /* Copy ati into atfill so we have locations of the next free space in atj */
2005: PetscMemcpy(atfill,ati,an*sizeof(PetscInt));
2007: /* Walk through A row-wise and mark nonzero entries of A^T. */
2008: for (i=0;i<am;i++) {
2009: anzj = ai[i+1] - ai[i];
2010: for (j=0;j<anzj;j++) {
2011: atj[atfill[*aj]] = i;
2012: atfill[*aj++] += 1;
2013: }
2014: }
2016: /* Clean up temporary space and complete requests. */
2017: PetscFree(atfill);
2018: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);
2019: MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
2021: b = (Mat_SeqAIJ*)((*B)->data);
2022: b->free_a = PETSC_FALSE;
2023: b->free_ij = PETSC_TRUE;
2024: b->nonew = 0;
2025: return(0);
2026: }
2030: PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B)
2031: {
2032: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2033: Mat C;
2035: PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col;
2036: MatScalar *array = a->a;
2039: if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place");
2041: if (reuse == MAT_INITIAL_MATRIX || *B == A) {
2042: PetscCalloc1(1+A->cmap->n,&col);
2044: for (i=0; i<ai[m]; i++) col[aj[i]] += 1;
2045: MatCreate(PetscObjectComm((PetscObject)A),&C);
2046: MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);
2047: MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
2048: MatSetType(C,((PetscObject)A)->type_name);
2049: MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);
2050: PetscFree(col);
2051: } else {
2052: C = *B;
2053: }
2055: for (i=0; i<m; i++) {
2056: len = ai[i+1]-ai[i];
2057: MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);
2058: array += len;
2059: aj += len;
2060: }
2061: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2062: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2064: if (reuse == MAT_INITIAL_MATRIX || *B != A) {
2065: *B = C;
2066: } else {
2067: MatHeaderMerge(A,C);
2068: }
2069: return(0);
2070: }
2074: PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2075: {
2076: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data;
2077: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2078: MatScalar *va,*vb;
2080: PetscInt ma,na,mb,nb, i;
2083: bij = (Mat_SeqAIJ*) B->data;
2085: MatGetSize(A,&ma,&na);
2086: MatGetSize(B,&mb,&nb);
2087: if (ma!=nb || na!=mb) {
2088: *f = PETSC_FALSE;
2089: return(0);
2090: }
2091: aii = aij->i; bii = bij->i;
2092: adx = aij->j; bdx = bij->j;
2093: va = aij->a; vb = bij->a;
2094: PetscMalloc1(ma,&aptr);
2095: PetscMalloc1(mb,&bptr);
2096: for (i=0; i<ma; i++) aptr[i] = aii[i];
2097: for (i=0; i<mb; i++) bptr[i] = bii[i];
2099: *f = PETSC_TRUE;
2100: for (i=0; i<ma; i++) {
2101: while (aptr[i]<aii[i+1]) {
2102: PetscInt idc,idr;
2103: PetscScalar vc,vr;
2104: /* column/row index/value */
2105: idc = adx[aptr[i]];
2106: idr = bdx[bptr[idc]];
2107: vc = va[aptr[i]];
2108: vr = vb[bptr[idc]];
2109: if (i!=idr || PetscAbsScalar(vc-vr) > tol) {
2110: *f = PETSC_FALSE;
2111: goto done;
2112: } else {
2113: aptr[i]++;
2114: if (B || i!=idc) bptr[idc]++;
2115: }
2116: }
2117: }
2118: done:
2119: PetscFree(aptr);
2120: PetscFree(bptr);
2121: return(0);
2122: }
2126: PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
2127: {
2128: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data;
2129: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
2130: MatScalar *va,*vb;
2132: PetscInt ma,na,mb,nb, i;
2135: bij = (Mat_SeqAIJ*) B->data;
2137: MatGetSize(A,&ma,&na);
2138: MatGetSize(B,&mb,&nb);
2139: if (ma!=nb || na!=mb) {
2140: *f = PETSC_FALSE;
2141: return(0);
2142: }
2143: aii = aij->i; bii = bij->i;
2144: adx = aij->j; bdx = bij->j;
2145: va = aij->a; vb = bij->a;
2146: PetscMalloc1(ma,&aptr);
2147: PetscMalloc1(mb,&bptr);
2148: for (i=0; i<ma; i++) aptr[i] = aii[i];
2149: for (i=0; i<mb; i++) bptr[i] = bii[i];
2151: *f = PETSC_TRUE;
2152: for (i=0; i<ma; i++) {
2153: while (aptr[i]<aii[i+1]) {
2154: PetscInt idc,idr;
2155: PetscScalar vc,vr;
2156: /* column/row index/value */
2157: idc = adx[aptr[i]];
2158: idr = bdx[bptr[idc]];
2159: vc = va[aptr[i]];
2160: vr = vb[bptr[idc]];
2161: if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) {
2162: *f = PETSC_FALSE;
2163: goto done;
2164: } else {
2165: aptr[i]++;
2166: if (B || i!=idc) bptr[idc]++;
2167: }
2168: }
2169: }
2170: done:
2171: PetscFree(aptr);
2172: PetscFree(bptr);
2173: return(0);
2174: }
2178: PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2179: {
2183: MatIsTranspose_SeqAIJ(A,A,tol,f);
2184: return(0);
2185: }
2189: PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2190: {
2194: MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);
2195: return(0);
2196: }
2200: PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr)
2201: {
2202: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2203: PetscScalar *l,*r,x;
2204: MatScalar *v;
2206: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj;
2209: if (ll) {
2210: /* The local size is used so that VecMPI can be passed to this routine
2211: by MatDiagonalScale_MPIAIJ */
2212: VecGetLocalSize(ll,&m);
2213: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length");
2214: VecGetArray(ll,&l);
2215: v = a->a;
2216: for (i=0; i<m; i++) {
2217: x = l[i];
2218: M = a->i[i+1] - a->i[i];
2219: for (j=0; j<M; j++) (*v++) *= x;
2220: }
2221: VecRestoreArray(ll,&l);
2222: PetscLogFlops(nz);
2223: }
2224: if (rr) {
2225: VecGetLocalSize(rr,&n);
2226: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length");
2227: VecGetArray(rr,&r);
2228: v = a->a; jj = a->j;
2229: for (i=0; i<nz; i++) (*v++) *= r[*jj++];
2230: VecRestoreArray(rr,&r);
2231: PetscLogFlops(nz);
2232: }
2233: MatSeqAIJInvalidateDiagonal(A);
2234: return(0);
2235: }
2239: PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B)
2240: {
2241: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c;
2243: PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens;
2244: PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi;
2245: const PetscInt *irow,*icol;
2246: PetscInt nrows,ncols;
2247: PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen;
2248: MatScalar *a_new,*mat_a;
2249: Mat C;
2250: PetscBool stride;
2254: ISGetIndices(isrow,&irow);
2255: ISGetLocalSize(isrow,&nrows);
2256: ISGetLocalSize(iscol,&ncols);
2258: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
2259: if (stride) {
2260: ISStrideGetInfo(iscol,&first,&step);
2261: } else {
2262: first = 0;
2263: step = 0;
2264: }
2265: if (stride && step == 1) {
2266: /* special case of contiguous rows */
2267: PetscMalloc2(nrows,&lens,nrows,&starts);
2268: /* loop over new rows determining lens and starting points */
2269: for (i=0; i<nrows; i++) {
2270: kstart = ai[irow[i]];
2271: kend = kstart + ailen[irow[i]];
2272: starts[i] = kstart;
2273: for (k=kstart; k<kend; k++) {
2274: if (aj[k] >= first) {
2275: starts[i] = k;
2276: break;
2277: }
2278: }
2279: sum = 0;
2280: while (k < kend) {
2281: if (aj[k++] >= first+ncols) break;
2282: sum++;
2283: }
2284: lens[i] = sum;
2285: }
2286: /* create submatrix */
2287: if (scall == MAT_REUSE_MATRIX) {
2288: PetscInt n_cols,n_rows;
2289: MatGetSize(*B,&n_rows,&n_cols);
2290: if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
2291: MatZeroEntries(*B);
2292: C = *B;
2293: } else {
2294: PetscInt rbs,cbs;
2295: MatCreate(PetscObjectComm((PetscObject)A),&C);
2296: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2297: ISGetBlockSize(isrow,&rbs);
2298: ISGetBlockSize(iscol,&cbs);
2299: MatSetBlockSizes(C,rbs,cbs);
2300: MatSetType(C,((PetscObject)A)->type_name);
2301: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2302: }
2303: c = (Mat_SeqAIJ*)C->data;
2305: /* loop over rows inserting into submatrix */
2306: a_new = c->a;
2307: j_new = c->j;
2308: i_new = c->i;
2310: for (i=0; i<nrows; i++) {
2311: ii = starts[i];
2312: lensi = lens[i];
2313: for (k=0; k<lensi; k++) {
2314: *j_new++ = aj[ii+k] - first;
2315: }
2316: PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));
2317: a_new += lensi;
2318: i_new[i+1] = i_new[i] + lensi;
2319: c->ilen[i] = lensi;
2320: }
2321: PetscFree2(lens,starts);
2322: } else {
2323: ISGetIndices(iscol,&icol);
2324: PetscCalloc1(oldcols,&smap);
2325: PetscMalloc1(1+nrows,&lens);
2326: for (i=0; i<ncols; i++) {
2327: #if defined(PETSC_USE_DEBUG)
2328: if (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);
2329: #endif
2330: smap[icol[i]] = i+1;
2331: }
2333: /* determine lens of each row */
2334: for (i=0; i<nrows; i++) {
2335: kstart = ai[irow[i]];
2336: kend = kstart + a->ilen[irow[i]];
2337: lens[i] = 0;
2338: for (k=kstart; k<kend; k++) {
2339: if (smap[aj[k]]) {
2340: lens[i]++;
2341: }
2342: }
2343: }
2344: /* Create and fill new matrix */
2345: if (scall == MAT_REUSE_MATRIX) {
2346: PetscBool equal;
2348: c = (Mat_SeqAIJ*)((*B)->data);
2349: if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size");
2350: PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);
2351: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros");
2352: PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));
2353: C = *B;
2354: } else {
2355: PetscInt rbs,cbs;
2356: MatCreate(PetscObjectComm((PetscObject)A),&C);
2357: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2358: ISGetBlockSize(isrow,&rbs);
2359: ISGetBlockSize(iscol,&cbs);
2360: MatSetBlockSizes(C,rbs,cbs);
2361: MatSetType(C,((PetscObject)A)->type_name);
2362: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2363: }
2364: c = (Mat_SeqAIJ*)(C->data);
2365: for (i=0; i<nrows; i++) {
2366: row = irow[i];
2367: kstart = ai[row];
2368: kend = kstart + a->ilen[row];
2369: mat_i = c->i[i];
2370: mat_j = c->j + mat_i;
2371: mat_a = c->a + mat_i;
2372: mat_ilen = c->ilen + i;
2373: for (k=kstart; k<kend; k++) {
2374: if ((tcol=smap[a->j[k]])) {
2375: *mat_j++ = tcol - 1;
2376: *mat_a++ = a->a[k];
2377: (*mat_ilen)++;
2379: }
2380: }
2381: }
2382: /* Free work space */
2383: ISRestoreIndices(iscol,&icol);
2384: PetscFree(smap);
2385: PetscFree(lens);
2386: /* sort */
2387: for (i = 0; i < nrows; i++) {
2388: PetscInt ilen;
2390: mat_i = c->i[i];
2391: mat_j = c->j + mat_i;
2392: mat_a = c->a + mat_i;
2393: ilen = c->ilen[i];
2394: PetscSortIntWithMatScalarArray(ilen,mat_j,mat_a);
2395: }
2396: }
2397: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2398: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2400: ISRestoreIndices(isrow,&irow);
2401: *B = C;
2402: return(0);
2403: }
2407: PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat)
2408: {
2410: Mat B;
2413: if (scall == MAT_INITIAL_MATRIX) {
2414: MatCreate(subComm,&B);
2415: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);
2416: MatSetBlockSizesFromMats(B,mat,mat);
2417: MatSetType(B,MATSEQAIJ);
2418: MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);
2419: *subMat = B;
2420: } else {
2421: MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);
2422: }
2423: return(0);
2424: }
2428: PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2429: {
2430: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2432: Mat outA;
2433: PetscBool row_identity,col_identity;
2436: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu");
2438: ISIdentity(row,&row_identity);
2439: ISIdentity(col,&col_identity);
2441: outA = inA;
2442: outA->factortype = MAT_FACTOR_LU;
2444: PetscObjectReference((PetscObject)row);
2445: ISDestroy(&a->row);
2447: a->row = row;
2449: PetscObjectReference((PetscObject)col);
2450: ISDestroy(&a->col);
2452: a->col = col;
2454: /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */
2455: ISDestroy(&a->icol);
2456: ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2457: PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);
2459: if (!a->solve_work) { /* this matrix may have been factored before */
2460: PetscMalloc1(inA->rmap->n+1,&a->solve_work);
2461: PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));
2462: }
2464: MatMarkDiagonal_SeqAIJ(inA);
2465: if (row_identity && col_identity) {
2466: MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);
2467: } else {
2468: MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);
2469: }
2470: return(0);
2471: }
2475: PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha)
2476: {
2477: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2478: PetscScalar oalpha = alpha;
2480: PetscBLASInt one = 1,bnz;
2483: PetscBLASIntCast(a->nz,&bnz);
2484: PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one));
2485: PetscLogFlops(a->nz);
2486: MatSeqAIJInvalidateDiagonal(inA);
2487: return(0);
2488: }
2492: PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2493: {
2495: PetscInt i;
2498: if (scall == MAT_INITIAL_MATRIX) {
2499: PetscMalloc1(n+1,B);
2500: }
2502: for (i=0; i<n; i++) {
2503: MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2504: }
2505: return(0);
2506: }
2510: PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov)
2511: {
2512: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2514: PetscInt row,i,j,k,l,m,n,*nidx,isz,val;
2515: const PetscInt *idx;
2516: PetscInt start,end,*ai,*aj;
2517: PetscBT table;
2520: m = A->rmap->n;
2521: ai = a->i;
2522: aj = a->j;
2524: if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used");
2526: PetscMalloc1(m+1,&nidx);
2527: PetscBTCreate(m,&table);
2529: for (i=0; i<is_max; i++) {
2530: /* Initialize the two local arrays */
2531: isz = 0;
2532: PetscBTMemzero(m,table);
2534: /* Extract the indices, assume there can be duplicate entries */
2535: ISGetIndices(is[i],&idx);
2536: ISGetLocalSize(is[i],&n);
2538: /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */
2539: for (j=0; j<n; ++j) {
2540: if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j];
2541: }
2542: ISRestoreIndices(is[i],&idx);
2543: ISDestroy(&is[i]);
2545: k = 0;
2546: for (j=0; j<ov; j++) { /* for each overlap */
2547: n = isz;
2548: for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */
2549: row = nidx[k];
2550: start = ai[row];
2551: end = ai[row+1];
2552: for (l = start; l<end; l++) {
2553: val = aj[l];
2554: if (!PetscBTLookupSet(table,val)) nidx[isz++] = val;
2555: }
2556: }
2557: }
2558: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));
2559: }
2560: PetscBTDestroy(&table);
2561: PetscFree(nidx);
2562: return(0);
2563: }
2565: /* -------------------------------------------------------------- */
2568: PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B)
2569: {
2570: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2572: PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n;
2573: const PetscInt *row,*col;
2574: PetscInt *cnew,j,*lens;
2575: IS icolp,irowp;
2576: PetscInt *cwork = NULL;
2577: PetscScalar *vwork = NULL;
2580: ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);
2581: ISGetIndices(irowp,&row);
2582: ISInvertPermutation(colp,PETSC_DECIDE,&icolp);
2583: ISGetIndices(icolp,&col);
2585: /* determine lengths of permuted rows */
2586: PetscMalloc1(m+1,&lens);
2587: for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i];
2588: MatCreate(PetscObjectComm((PetscObject)A),B);
2589: MatSetSizes(*B,m,n,m,n);
2590: MatSetBlockSizesFromMats(*B,A,A);
2591: MatSetType(*B,((PetscObject)A)->type_name);
2592: MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);
2593: PetscFree(lens);
2595: PetscMalloc1(n,&cnew);
2596: for (i=0; i<m; i++) {
2597: MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2598: for (j=0; j<nz; j++) cnew[j] = col[cwork[j]];
2599: MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);
2600: MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2601: }
2602: PetscFree(cnew);
2604: (*B)->assembled = PETSC_FALSE;
2606: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
2607: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
2608: ISRestoreIndices(irowp,&row);
2609: ISRestoreIndices(icolp,&col);
2610: ISDestroy(&irowp);
2611: ISDestroy(&icolp);
2612: return(0);
2613: }
2617: PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str)
2618: {
2622: /* If the two matrices have the same copy implementation, use fast copy. */
2623: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2624: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2625: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
2627: if (a->i[A->rmap->n] != b->i[B->rmap->n]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different");
2628: PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));
2629: } else {
2630: MatCopy_Basic(A,B,str);
2631: }
2632: return(0);
2633: }
2637: PetscErrorCode MatSetUp_SeqAIJ(Mat A)
2638: {
2642: MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);
2643: return(0);
2644: }
2648: PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[])
2649: {
2650: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2653: *array = a->a;
2654: return(0);
2655: }
2659: PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[])
2660: {
2662: return(0);
2663: }
2665: /*
2666: Computes the number of nonzeros per row needed for preallocation when X and Y
2667: have different nonzero structure.
2668: */
2671: PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz)
2672: {
2673: PetscInt i,j,k,nzx,nzy;
2676: /* Set the number of nonzeros in the new matrix */
2677: for (i=0; i<m; i++) {
2678: const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i];
2679: nzx = xi[i+1] - xi[i];
2680: nzy = yi[i+1] - yi[i];
2681: nnz[i] = 0;
2682: for (j=0,k=0; j<nzx; j++) { /* Point in X */
2683: for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */
2684: if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */
2685: nnz[i]++;
2686: }
2687: for (; k<nzy; k++) nnz[i]++;
2688: }
2689: return(0);
2690: }
2694: PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz)
2695: {
2696: PetscInt m = Y->rmap->N;
2697: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data;
2698: Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data;
2702: /* Set the number of nonzeros in the new matrix */
2703: MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);
2704: return(0);
2705: }
2709: PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
2710: {
2712: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data;
2713: PetscBLASInt one=1,bnz;
2716: PetscBLASIntCast(x->nz,&bnz);
2717: if (str == SAME_NONZERO_PATTERN) {
2718: PetscScalar alpha = a;
2719: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
2720: MatSeqAIJInvalidateDiagonal(Y);
2721: PetscObjectStateIncrease((PetscObject)Y);
2722: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
2723: MatAXPY_Basic(Y,a,X,str);
2724: } else {
2725: Mat B;
2726: PetscInt *nnz;
2727: PetscMalloc1(Y->rmap->N,&nnz);
2728: MatCreate(PetscObjectComm((PetscObject)Y),&B);
2729: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
2730: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
2731: MatSetBlockSizesFromMats(B,Y,Y);
2732: MatSetType(B,(MatType) ((PetscObject)Y)->type_name);
2733: MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);
2734: MatSeqAIJSetPreallocation(B,0,nnz);
2735: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2736: MatHeaderReplace(Y,B);
2737: PetscFree(nnz);
2738: }
2739: return(0);
2740: }
2744: PetscErrorCode MatConjugate_SeqAIJ(Mat mat)
2745: {
2746: #if defined(PETSC_USE_COMPLEX)
2747: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
2748: PetscInt i,nz;
2749: PetscScalar *a;
2752: nz = aij->nz;
2753: a = aij->a;
2754: for (i=0; i<nz; i++) a[i] = PetscConj(a[i]);
2755: #else
2757: #endif
2758: return(0);
2759: }
2763: PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2764: {
2765: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2767: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2768: PetscReal atmp;
2769: PetscScalar *x;
2770: MatScalar *aa;
2773: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2774: aa = a->a;
2775: ai = a->i;
2776: aj = a->j;
2778: VecSet(v,0.0);
2779: VecGetArray(v,&x);
2780: VecGetLocalSize(v,&n);
2781: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2782: for (i=0; i<m; i++) {
2783: ncols = ai[1] - ai[0]; ai++;
2784: x[i] = 0.0;
2785: for (j=0; j<ncols; j++) {
2786: atmp = PetscAbsScalar(*aa);
2787: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
2788: aa++; aj++;
2789: }
2790: }
2791: VecRestoreArray(v,&x);
2792: return(0);
2793: }
2797: PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2798: {
2799: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2801: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2802: PetscScalar *x;
2803: MatScalar *aa;
2806: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2807: aa = a->a;
2808: ai = a->i;
2809: aj = a->j;
2811: VecSet(v,0.0);
2812: VecGetArray(v,&x);
2813: VecGetLocalSize(v,&n);
2814: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2815: for (i=0; i<m; i++) {
2816: ncols = ai[1] - ai[0]; ai++;
2817: if (ncols == A->cmap->n) { /* row is dense */
2818: x[i] = *aa; if (idx) idx[i] = 0;
2819: } else { /* row is sparse so already KNOW maximum is 0.0 or higher */
2820: x[i] = 0.0;
2821: if (idx) {
2822: idx[i] = 0; /* in case ncols is zero */
2823: for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */
2824: if (aj[j] > j) {
2825: idx[i] = j;
2826: break;
2827: }
2828: }
2829: }
2830: }
2831: for (j=0; j<ncols; j++) {
2832: if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
2833: aa++; aj++;
2834: }
2835: }
2836: VecRestoreArray(v,&x);
2837: return(0);
2838: }
2842: PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2843: {
2844: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2846: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2847: PetscReal atmp;
2848: PetscScalar *x;
2849: MatScalar *aa;
2852: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2853: aa = a->a;
2854: ai = a->i;
2855: aj = a->j;
2857: VecSet(v,0.0);
2858: VecGetArray(v,&x);
2859: VecGetLocalSize(v,&n);
2860: if (n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %D vs. %D rows", A->rmap->n, n);
2861: for (i=0; i<m; i++) {
2862: ncols = ai[1] - ai[0]; ai++;
2863: if (ncols) {
2864: /* Get first nonzero */
2865: for (j = 0; j < ncols; j++) {
2866: atmp = PetscAbsScalar(aa[j]);
2867: if (atmp > 1.0e-12) {
2868: x[i] = atmp;
2869: if (idx) idx[i] = aj[j];
2870: break;
2871: }
2872: }
2873: if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;}
2874: } else {
2875: x[i] = 0.0; if (idx) idx[i] = 0;
2876: }
2877: for (j = 0; j < ncols; j++) {
2878: atmp = PetscAbsScalar(*aa);
2879: if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
2880: aa++; aj++;
2881: }
2882: }
2883: VecRestoreArray(v,&x);
2884: return(0);
2885: }
2889: PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2890: {
2891: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2892: PetscErrorCode ierr;
2893: PetscInt i,j,m = A->rmap->n,ncols,n;
2894: const PetscInt *ai,*aj;
2895: PetscScalar *x;
2896: const MatScalar *aa;
2899: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2900: aa = a->a;
2901: ai = a->i;
2902: aj = a->j;
2904: VecSet(v,0.0);
2905: VecGetArray(v,&x);
2906: VecGetLocalSize(v,&n);
2907: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2908: for (i=0; i<m; i++) {
2909: ncols = ai[1] - ai[0]; ai++;
2910: if (ncols == A->cmap->n) { /* row is dense */
2911: x[i] = *aa; if (idx) idx[i] = 0;
2912: } else { /* row is sparse so already KNOW minimum is 0.0 or lower */
2913: x[i] = 0.0;
2914: if (idx) { /* find first implicit 0.0 in the row */
2915: idx[i] = 0; /* in case ncols is zero */
2916: for (j=0; j<ncols; j++) {
2917: if (aj[j] > j) {
2918: idx[i] = j;
2919: break;
2920: }
2921: }
2922: }
2923: }
2924: for (j=0; j<ncols; j++) {
2925: if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
2926: aa++; aj++;
2927: }
2928: }
2929: VecRestoreArray(v,&x);
2930: return(0);
2931: }
2933: #include <petscblaslapack.h>
2934: #include <petsc/private/kernels/blockinvert.h>
2938: PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values)
2939: {
2940: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
2942: PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j;
2943: MatScalar *diag,work[25],*v_work;
2944: PetscReal shift = 0.0;
2947: if (a->ibdiagvalid) {
2948: if (values) *values = a->ibdiag;
2949: return(0);
2950: }
2951: MatMarkDiagonal_SeqAIJ(A);
2952: if (!a->ibdiag) {
2953: PetscMalloc1(bs2*mbs,&a->ibdiag);
2954: PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));
2955: }
2956: diag = a->ibdiag;
2957: if (values) *values = a->ibdiag;
2958: /* factor and invert each block */
2959: switch (bs) {
2960: case 1:
2961: for (i=0; i<mbs; i++) {
2962: MatGetValues(A,1,&i,1,&i,diag+i);
2963: diag[i] = (PetscScalar)1.0 / (diag[i] + shift);
2964: }
2965: break;
2966: case 2:
2967: for (i=0; i<mbs; i++) {
2968: ij[0] = 2*i; ij[1] = 2*i + 1;
2969: MatGetValues(A,2,ij,2,ij,diag);
2970: PetscKernel_A_gets_inverse_A_2(diag,shift);
2971: PetscKernel_A_gets_transpose_A_2(diag);
2972: diag += 4;
2973: }
2974: break;
2975: case 3:
2976: for (i=0; i<mbs; i++) {
2977: ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2;
2978: MatGetValues(A,3,ij,3,ij,diag);
2979: PetscKernel_A_gets_inverse_A_3(diag,shift);
2980: PetscKernel_A_gets_transpose_A_3(diag);
2981: diag += 9;
2982: }
2983: break;
2984: case 4:
2985: for (i=0; i<mbs; i++) {
2986: ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3;
2987: MatGetValues(A,4,ij,4,ij,diag);
2988: PetscKernel_A_gets_inverse_A_4(diag,shift);
2989: PetscKernel_A_gets_transpose_A_4(diag);
2990: diag += 16;
2991: }
2992: break;
2993: case 5:
2994: for (i=0; i<mbs; i++) {
2995: ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4;
2996: MatGetValues(A,5,ij,5,ij,diag);
2997: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);
2998: PetscKernel_A_gets_transpose_A_5(diag);
2999: diag += 25;
3000: }
3001: break;
3002: case 6:
3003: for (i=0; i<mbs; i++) {
3004: 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;
3005: MatGetValues(A,6,ij,6,ij,diag);
3006: PetscKernel_A_gets_inverse_A_6(diag,shift);
3007: PetscKernel_A_gets_transpose_A_6(diag);
3008: diag += 36;
3009: }
3010: break;
3011: case 7:
3012: for (i=0; i<mbs; i++) {
3013: 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;
3014: MatGetValues(A,7,ij,7,ij,diag);
3015: PetscKernel_A_gets_inverse_A_7(diag,shift);
3016: PetscKernel_A_gets_transpose_A_7(diag);
3017: diag += 49;
3018: }
3019: break;
3020: default:
3021: PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);
3022: for (i=0; i<mbs; i++) {
3023: for (j=0; j<bs; j++) {
3024: IJ[j] = bs*i + j;
3025: }
3026: MatGetValues(A,bs,IJ,bs,IJ,diag);
3027: PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);
3028: PetscKernel_A_gets_transpose_A_N(diag,bs);
3029: diag += bs2;
3030: }
3031: PetscFree3(v_work,v_pivots,IJ);
3032: }
3033: a->ibdiagvalid = PETSC_TRUE;
3034: return(0);
3035: }
3039: static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx)
3040: {
3042: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3043: PetscScalar a;
3044: PetscInt m,n,i,j,col;
3047: if (!x->assembled) {
3048: MatGetSize(x,&m,&n);
3049: for (i=0; i<m; i++) {
3050: for (j=0; j<aij->imax[i]; j++) {
3051: PetscRandomGetValue(rctx,&a);
3052: col = (PetscInt)(n*PetscRealPart(a));
3053: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3054: }
3055: }
3056: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded");
3057: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3058: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3059: return(0);
3060: }
3064: PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a)
3065: {
3067: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data;
3070: if (!Y->preallocated || !aij->nz) {
3071: MatSeqAIJSetPreallocation(Y,1,NULL);
3072: }
3073: MatShift_Basic(Y,a);
3074: return(0);
3075: }
3077: /* -------------------------------------------------------------------*/
3078: static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ,
3079: MatGetRow_SeqAIJ,
3080: MatRestoreRow_SeqAIJ,
3081: MatMult_SeqAIJ,
3082: /* 4*/ MatMultAdd_SeqAIJ,
3083: MatMultTranspose_SeqAIJ,
3084: MatMultTransposeAdd_SeqAIJ,
3085: 0,
3086: 0,
3087: 0,
3088: /* 10*/ 0,
3089: MatLUFactor_SeqAIJ,
3090: 0,
3091: MatSOR_SeqAIJ,
3092: MatTranspose_SeqAIJ,
3093: /*1 5*/ MatGetInfo_SeqAIJ,
3094: MatEqual_SeqAIJ,
3095: MatGetDiagonal_SeqAIJ,
3096: MatDiagonalScale_SeqAIJ,
3097: MatNorm_SeqAIJ,
3098: /* 20*/ 0,
3099: MatAssemblyEnd_SeqAIJ,
3100: MatSetOption_SeqAIJ,
3101: MatZeroEntries_SeqAIJ,
3102: /* 24*/ MatZeroRows_SeqAIJ,
3103: 0,
3104: 0,
3105: 0,
3106: 0,
3107: /* 29*/ MatSetUp_SeqAIJ,
3108: 0,
3109: 0,
3110: 0,
3111: 0,
3112: /* 34*/ MatDuplicate_SeqAIJ,
3113: 0,
3114: 0,
3115: MatILUFactor_SeqAIJ,
3116: 0,
3117: /* 39*/ MatAXPY_SeqAIJ,
3118: MatGetSubMatrices_SeqAIJ,
3119: MatIncreaseOverlap_SeqAIJ,
3120: MatGetValues_SeqAIJ,
3121: MatCopy_SeqAIJ,
3122: /* 44*/ MatGetRowMax_SeqAIJ,
3123: MatScale_SeqAIJ,
3124: MatShift_SeqAIJ,
3125: MatDiagonalSet_SeqAIJ,
3126: MatZeroRowsColumns_SeqAIJ,
3127: /* 49*/ MatSetRandom_SeqAIJ,
3128: MatGetRowIJ_SeqAIJ,
3129: MatRestoreRowIJ_SeqAIJ,
3130: MatGetColumnIJ_SeqAIJ,
3131: MatRestoreColumnIJ_SeqAIJ,
3132: /* 54*/ MatFDColoringCreate_SeqXAIJ,
3133: 0,
3134: 0,
3135: MatPermute_SeqAIJ,
3136: 0,
3137: /* 59*/ 0,
3138: MatDestroy_SeqAIJ,
3139: MatView_SeqAIJ,
3140: 0,
3141: MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ,
3142: /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ,
3143: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ,
3144: 0,
3145: 0,
3146: 0,
3147: /* 69*/ MatGetRowMaxAbs_SeqAIJ,
3148: MatGetRowMinAbs_SeqAIJ,
3149: 0,
3150: MatSetColoring_SeqAIJ,
3151: 0,
3152: /* 74*/ MatSetValuesAdifor_SeqAIJ,
3153: MatFDColoringApply_AIJ,
3154: 0,
3155: 0,
3156: 0,
3157: /* 79*/ MatFindZeroDiagonals_SeqAIJ,
3158: 0,
3159: 0,
3160: 0,
3161: MatLoad_SeqAIJ,
3162: /* 84*/ MatIsSymmetric_SeqAIJ,
3163: MatIsHermitian_SeqAIJ,
3164: 0,
3165: 0,
3166: 0,
3167: /* 89*/ MatMatMult_SeqAIJ_SeqAIJ,
3168: MatMatMultSymbolic_SeqAIJ_SeqAIJ,
3169: MatMatMultNumeric_SeqAIJ_SeqAIJ,
3170: MatPtAP_SeqAIJ_SeqAIJ,
3171: MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy,
3172: /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ,
3173: MatMatTransposeMult_SeqAIJ_SeqAIJ,
3174: MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ,
3175: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ,
3176: 0,
3177: /* 99*/ 0,
3178: 0,
3179: 0,
3180: MatConjugate_SeqAIJ,
3181: 0,
3182: /*104*/ MatSetValuesRow_SeqAIJ,
3183: MatRealPart_SeqAIJ,
3184: MatImaginaryPart_SeqAIJ,
3185: 0,
3186: 0,
3187: /*109*/ MatMatSolve_SeqAIJ,
3188: 0,
3189: MatGetRowMin_SeqAIJ,
3190: 0,
3191: MatMissingDiagonal_SeqAIJ,
3192: /*114*/ 0,
3193: 0,
3194: 0,
3195: 0,
3196: 0,
3197: /*119*/ 0,
3198: 0,
3199: 0,
3200: 0,
3201: MatGetMultiProcBlock_SeqAIJ,
3202: /*124*/ MatFindNonzeroRows_SeqAIJ,
3203: MatGetColumnNorms_SeqAIJ,
3204: MatInvertBlockDiagonal_SeqAIJ,
3205: 0,
3206: 0,
3207: /*129*/ 0,
3208: MatTransposeMatMult_SeqAIJ_SeqAIJ,
3209: MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ,
3210: MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ,
3211: MatTransposeColoringCreate_SeqAIJ,
3212: /*134*/ MatTransColoringApplySpToDen_SeqAIJ,
3213: MatTransColoringApplyDenToSp_SeqAIJ,
3214: MatRARt_SeqAIJ_SeqAIJ,
3215: MatRARtSymbolic_SeqAIJ_SeqAIJ,
3216: MatRARtNumeric_SeqAIJ_SeqAIJ,
3217: /*139*/0,
3218: 0,
3219: 0,
3220: MatFDColoringSetUp_SeqXAIJ,
3221: MatFindOffBlockDiagonalEntries_SeqAIJ,
3222: /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ
3223: };
3227: PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices)
3228: {
3229: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3230: PetscInt i,nz,n;
3233: nz = aij->maxnz;
3234: n = mat->rmap->n;
3235: for (i=0; i<nz; i++) {
3236: aij->j[i] = indices[i];
3237: }
3238: aij->nz = nz;
3239: for (i=0; i<n; i++) {
3240: aij->ilen[i] = aij->imax[i];
3241: }
3242: return(0);
3243: }
3247: /*@
3248: MatSeqAIJSetColumnIndices - Set the column indices for all the rows
3249: in the matrix.
3251: Input Parameters:
3252: + mat - the SeqAIJ matrix
3253: - indices - the column indices
3255: Level: advanced
3257: Notes:
3258: This can be called if you have precomputed the nonzero structure of the
3259: matrix and want to provide it to the matrix object to improve the performance
3260: of the MatSetValues() operation.
3262: You MUST have set the correct numbers of nonzeros per row in the call to
3263: MatCreateSeqAIJ(), and the columns indices MUST be sorted.
3265: MUST be called before any calls to MatSetValues();
3267: The indices should start with zero, not one.
3269: @*/
3270: PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices)
3271: {
3277: PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
3278: return(0);
3279: }
3281: /* ----------------------------------------------------------------------------------------*/
3285: PetscErrorCode MatStoreValues_SeqAIJ(Mat mat)
3286: {
3287: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3289: size_t nz = aij->i[mat->rmap->n];
3292: if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3294: /* allocate space for values if not already there */
3295: if (!aij->saved_values) {
3296: PetscMalloc1(nz+1,&aij->saved_values);
3297: PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
3298: }
3300: /* copy values over */
3301: PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));
3302: return(0);
3303: }
3307: /*@
3308: MatStoreValues - Stashes a copy of the matrix values; this allows, for
3309: example, reuse of the linear part of a Jacobian, while recomputing the
3310: nonlinear portion.
3312: Collect on Mat
3314: Input Parameters:
3315: . mat - the matrix (currently only AIJ matrices support this option)
3317: Level: advanced
3319: Common Usage, with SNESSolve():
3320: $ Create Jacobian matrix
3321: $ Set linear terms into matrix
3322: $ Apply boundary conditions to matrix, at this time matrix must have
3323: $ final nonzero structure (i.e. setting the nonlinear terms and applying
3324: $ boundary conditions again will not change the nonzero structure
3325: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3326: $ MatStoreValues(mat);
3327: $ Call SNESSetJacobian() with matrix
3328: $ In your Jacobian routine
3329: $ MatRetrieveValues(mat);
3330: $ Set nonlinear terms in matrix
3332: Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself:
3333: $ // build linear portion of Jacobian
3334: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3335: $ MatStoreValues(mat);
3336: $ loop over nonlinear iterations
3337: $ MatRetrieveValues(mat);
3338: $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian
3339: $ // call MatAssemblyBegin/End() on matrix
3340: $ Solve linear system with Jacobian
3341: $ endloop
3343: Notes:
3344: Matrix must already be assemblied before calling this routine
3345: Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before
3346: calling this routine.
3348: When this is called multiple times it overwrites the previous set of stored values
3349: and does not allocated additional space.
3351: .seealso: MatRetrieveValues()
3353: @*/
3354: PetscErrorCode MatStoreValues(Mat mat)
3355: {
3360: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3361: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3362: PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));
3363: return(0);
3364: }
3368: PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat)
3369: {
3370: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3372: PetscInt nz = aij->i[mat->rmap->n];
3375: if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3376: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
3377: /* copy values over */
3378: PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));
3379: return(0);
3380: }
3384: /*@
3385: MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for
3386: example, reuse of the linear part of a Jacobian, while recomputing the
3387: nonlinear portion.
3389: Collect on Mat
3391: Input Parameters:
3392: . mat - the matrix (currently on AIJ matrices support this option)
3394: Level: advanced
3396: .seealso: MatStoreValues()
3398: @*/
3399: PetscErrorCode MatRetrieveValues(Mat mat)
3400: {
3405: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3406: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3407: PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));
3408: return(0);
3409: }
3412: /* --------------------------------------------------------------------------------*/
3415: /*@C
3416: MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format
3417: (the default parallel PETSc format). For good matrix assembly performance
3418: the user should preallocate the matrix storage by setting the parameter nz
3419: (or the array nnz). By setting these parameters accurately, performance
3420: during matrix assembly can be increased by more than a factor of 50.
3422: Collective on MPI_Comm
3424: Input Parameters:
3425: + comm - MPI communicator, set to PETSC_COMM_SELF
3426: . m - number of rows
3427: . n - number of columns
3428: . nz - number of nonzeros per row (same for all rows)
3429: - nnz - array containing the number of nonzeros in the various rows
3430: (possibly different for each row) or NULL
3432: Output Parameter:
3433: . A - the matrix
3435: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3436: MatXXXXSetPreallocation() paradgm instead of this routine directly.
3437: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
3439: Notes:
3440: If nnz is given then nz is ignored
3442: The AIJ format (also called the Yale sparse matrix format or
3443: compressed row storage), is fully compatible with standard Fortran 77
3444: storage. That is, the stored row and column indices can begin at
3445: either one (as in Fortran) or zero. See the users' manual for details.
3447: Specify the preallocated storage with either nz or nnz (not both).
3448: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3449: allocation. For large problems you MUST preallocate memory or you
3450: will get TERRIBLE performance, see the users' manual chapter on matrices.
3452: By default, this format uses inodes (identical nodes) when possible, to
3453: improve numerical efficiency of matrix-vector products and solves. We
3454: search for consecutive rows with the same nonzero structure, thereby
3455: reusing matrix information to achieve increased efficiency.
3457: Options Database Keys:
3458: + -mat_no_inode - Do not use inodes
3459: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3461: Level: intermediate
3463: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays()
3465: @*/
3466: PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3467: {
3471: MatCreate(comm,A);
3472: MatSetSizes(*A,m,n,m,n);
3473: MatSetType(*A,MATSEQAIJ);
3474: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
3475: return(0);
3476: }
3480: /*@C
3481: MatSeqAIJSetPreallocation - For good matrix assembly performance
3482: the user should preallocate the matrix storage by setting the parameter nz
3483: (or the array nnz). By setting these parameters accurately, performance
3484: during matrix assembly can be increased by more than a factor of 50.
3486: Collective on MPI_Comm
3488: Input Parameters:
3489: + B - The matrix
3490: . nz - number of nonzeros per row (same for all rows)
3491: - nnz - array containing the number of nonzeros in the various rows
3492: (possibly different for each row) or NULL
3494: Notes:
3495: If nnz is given then nz is ignored
3497: The AIJ format (also called the Yale sparse matrix format or
3498: compressed row storage), is fully compatible with standard Fortran 77
3499: storage. That is, the stored row and column indices can begin at
3500: either one (as in Fortran) or zero. See the users' manual for details.
3502: Specify the preallocated storage with either nz or nnz (not both).
3503: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3504: allocation. For large problems you MUST preallocate memory or you
3505: will get TERRIBLE performance, see the users' manual chapter on matrices.
3507: You can call MatGetInfo() to get information on how effective the preallocation was;
3508: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3509: You can also run with the option -info and look for messages with the string
3510: malloc in them to see if additional memory allocation was needed.
3512: Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix
3513: entries or columns indices
3515: By default, this format uses inodes (identical nodes) when possible, to
3516: improve numerical efficiency of matrix-vector products and solves. We
3517: search for consecutive rows with the same nonzero structure, thereby
3518: reusing matrix information to achieve increased efficiency.
3520: Options Database Keys:
3521: + -mat_no_inode - Do not use inodes
3522: . -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3523: - -mat_aij_oneindex - Internally use indexing starting at 1
3524: rather than 0. Note that when calling MatSetValues(),
3525: the user still MUST index entries starting at 0!
3527: Level: intermediate
3529: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo()
3531: @*/
3532: PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[])
3533: {
3539: PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));
3540: return(0);
3541: }
3545: PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz)
3546: {
3547: Mat_SeqAIJ *b;
3548: PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;
3550: PetscInt i;
3553: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3554: if (nz == MAT_SKIP_ALLOCATION) {
3555: skipallocation = PETSC_TRUE;
3556: nz = 0;
3557: }
3559: PetscLayoutSetUp(B->rmap);
3560: PetscLayoutSetUp(B->cmap);
3562: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3563: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
3564: if (nnz) {
3565: for (i=0; i<B->rmap->n; i++) {
3566: 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]);
3567: 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);
3568: }
3569: }
3571: B->preallocated = PETSC_TRUE;
3573: b = (Mat_SeqAIJ*)B->data;
3575: if (!skipallocation) {
3576: if (!b->imax) {
3577: PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);
3578: PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));
3579: }
3580: if (!nnz) {
3581: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10;
3582: else if (nz < 0) nz = 1;
3583: for (i=0; i<B->rmap->n; i++) b->imax[i] = nz;
3584: nz = nz*B->rmap->n;
3585: } else {
3586: nz = 0;
3587: for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];}
3588: }
3589: /* b->ilen will count nonzeros in each row so far. */
3590: for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0;
3592: /* allocate the matrix space */
3593: /* FIXME: should B's old memory be unlogged? */
3594: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
3595: PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);
3596: PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
3597: b->i[0] = 0;
3598: for (i=1; i<B->rmap->n+1; i++) {
3599: b->i[i] = b->i[i-1] + b->imax[i-1];
3600: }
3601: b->singlemalloc = PETSC_TRUE;
3602: b->free_a = PETSC_TRUE;
3603: b->free_ij = PETSC_TRUE;
3604: } else {
3605: b->free_a = PETSC_FALSE;
3606: b->free_ij = PETSC_FALSE;
3607: }
3609: b->nz = 0;
3610: b->maxnz = nz;
3611: B->info.nz_unneeded = (double)b->maxnz;
3612: if (realalloc) {
3613: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
3614: }
3615: return(0);
3616: }
3618: #undef __FUNCT__
3620: /*@
3621: MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format.
3623: Input Parameters:
3624: + B - the matrix
3625: . i - the indices into j for the start of each row (starts with zero)
3626: . j - the column indices for each row (starts with zero) these must be sorted for each row
3627: - v - optional values in the matrix
3629: Level: developer
3631: The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays()
3633: .keywords: matrix, aij, compressed row, sparse, sequential
3635: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ
3636: @*/
3637: PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[])
3638: {
3644: PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));
3645: return(0);
3646: }
3648: #undef __FUNCT__
3650: PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[])
3651: {
3652: PetscInt i;
3653: PetscInt m,n;
3654: PetscInt nz;
3655: PetscInt *nnz, nz_max = 0;
3656: PetscScalar *values;
3660: if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]);
3662: PetscLayoutSetUp(B->rmap);
3663: PetscLayoutSetUp(B->cmap);
3665: MatGetSize(B, &m, &n);
3666: PetscMalloc1(m+1, &nnz);
3667: for (i = 0; i < m; i++) {
3668: nz = Ii[i+1]- Ii[i];
3669: nz_max = PetscMax(nz_max, nz);
3670: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz);
3671: nnz[i] = nz;
3672: }
3673: MatSeqAIJSetPreallocation(B, 0, nnz);
3674: PetscFree(nnz);
3676: if (v) {
3677: values = (PetscScalar*) v;
3678: } else {
3679: PetscCalloc1(nz_max, &values);
3680: }
3682: for (i = 0; i < m; i++) {
3683: nz = Ii[i+1] - Ii[i];
3684: MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);
3685: }
3687: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3688: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3690: if (!v) {
3691: PetscFree(values);
3692: }
3693: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3694: return(0);
3695: }
3697: #include <../src/mat/impls/dense/seq/dense.h>
3698: #include <petsc/private/kernels/petscaxpy.h>
3702: /*
3703: Computes (B'*A')' since computing B*A directly is untenable
3705: n p p
3706: ( ) ( ) ( )
3707: m ( A ) * n ( B ) = m ( C )
3708: ( ) ( ) ( )
3710: */
3711: PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C)
3712: {
3713: PetscErrorCode ierr;
3714: Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data;
3715: Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data;
3716: Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data;
3717: PetscInt i,n,m,q,p;
3718: const PetscInt *ii,*idx;
3719: const PetscScalar *b,*a,*a_q;
3720: PetscScalar *c,*c_q;
3723: m = A->rmap->n;
3724: n = A->cmap->n;
3725: p = B->cmap->n;
3726: a = sub_a->v;
3727: b = sub_b->a;
3728: c = sub_c->v;
3729: PetscMemzero(c,m*p*sizeof(PetscScalar));
3731: ii = sub_b->i;
3732: idx = sub_b->j;
3733: for (i=0; i<n; i++) {
3734: q = ii[i+1] - ii[i];
3735: while (q-->0) {
3736: c_q = c + m*(*idx);
3737: a_q = a + m*i;
3738: PetscKernelAXPY(c_q,*b,a_q,m);
3739: idx++;
3740: b++;
3741: }
3742: }
3743: return(0);
3744: }
3748: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
3749: {
3751: PetscInt m=A->rmap->n,n=B->cmap->n;
3752: Mat Cmat;
3755: 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);
3756: MatCreate(PetscObjectComm((PetscObject)A),&Cmat);
3757: MatSetSizes(Cmat,m,n,m,n);
3758: MatSetBlockSizesFromMats(Cmat,A,B);
3759: MatSetType(Cmat,MATSEQDENSE);
3760: MatSeqDenseSetPreallocation(Cmat,NULL);
3762: Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ;
3764: *C = Cmat;
3765: return(0);
3766: }
3768: /* ----------------------------------------------------------------*/
3771: PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
3772: {
3776: if (scall == MAT_INITIAL_MATRIX) {
3777: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
3778: MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);
3779: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
3780: }
3781: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
3782: MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);
3783: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
3784: return(0);
3785: }
3788: /*MC
3789: MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices,
3790: based on compressed sparse row format.
3792: Options Database Keys:
3793: . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions()
3795: Level: beginner
3797: .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType
3798: M*/
3800: /*MC
3801: MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices.
3803: This matrix type is identical to MATSEQAIJ when constructed with a single process communicator,
3804: and MATMPIAIJ otherwise. As a result, for single process communicators,
3805: MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported
3806: for communicators controlling multiple processes. It is recommended that you call both of
3807: the above preallocation routines for simplicity.
3809: Options Database Keys:
3810: . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions()
3812: Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when
3813: enough exist.
3815: Level: beginner
3817: .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ
3818: M*/
3820: /*MC
3821: MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices.
3823: This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator,
3824: and MATMPIAIJCRL otherwise. As a result, for single process communicators,
3825: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
3826: for communicators controlling multiple processes. It is recommended that you call both of
3827: the above preallocation routines for simplicity.
3829: Options Database Keys:
3830: . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions()
3832: Level: beginner
3834: .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL
3835: M*/
3837: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*);
3838: #if defined(PETSC_HAVE_ELEMENTAL)
3839: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
3840: #endif
3841: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*);
3843: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3844: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*);
3845: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*);
3846: #endif
3851: /*@C
3852: MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored
3854: Not Collective
3856: Input Parameter:
3857: . mat - a MATSEQAIJ matrix
3859: Output Parameter:
3860: . array - pointer to the data
3862: Level: intermediate
3864: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
3865: @*/
3866: PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array)
3867: {
3871: PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));
3872: return(0);
3873: }
3877: /*@C
3878: MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row
3880: Not Collective
3882: Input Parameter:
3883: . mat - a MATSEQAIJ matrix
3885: Output Parameter:
3886: . nz - the maximum number of nonzeros in any row
3888: Level: intermediate
3890: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
3891: @*/
3892: PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz)
3893: {
3894: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
3897: *nz = aij->rmax;
3898: return(0);
3899: }
3903: /*@C
3904: MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray()
3906: Not Collective
3908: Input Parameters:
3909: . mat - a MATSEQAIJ matrix
3910: . array - pointer to the data
3912: Level: intermediate
3914: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90()
3915: @*/
3916: PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array)
3917: {
3921: PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
3922: return(0);
3923: }
3927: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B)
3928: {
3929: Mat_SeqAIJ *b;
3931: PetscMPIInt size;
3934: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3935: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1");
3937: PetscNewLog(B,&b);
3939: B->data = (void*)b;
3941: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3943: b->row = 0;
3944: b->col = 0;
3945: b->icol = 0;
3946: b->reallocs = 0;
3947: b->ignorezeroentries = PETSC_FALSE;
3948: b->roworiented = PETSC_TRUE;
3949: b->nonew = 0;
3950: b->diag = 0;
3951: b->solve_work = 0;
3952: B->spptr = 0;
3953: b->saved_values = 0;
3954: b->idiag = 0;
3955: b->mdiag = 0;
3956: b->ssor_work = 0;
3957: b->omega = 1.0;
3958: b->fshift = 0.0;
3959: b->idiagvalid = PETSC_FALSE;
3960: b->ibdiagvalid = PETSC_FALSE;
3961: b->keepnonzeropattern = PETSC_FALSE;
3963: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
3964: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);
3965: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);
3967: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3968: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);
3969: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);
3970: #endif
3972: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);
3973: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);
3974: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);
3975: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);
3976: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);
3977: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);
3978: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);
3979: #if defined(PETSC_HAVE_ELEMENTAL)
3980: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);
3981: #endif
3982: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);
3983: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);
3984: PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);
3985: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);
3986: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);
3987: PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);
3988: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);
3989: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);
3990: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);
3991: MatCreate_SeqAIJ_Inode(B);
3992: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
3993: return(0);
3994: }
3998: /*
3999: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
4000: */
4001: PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
4002: {
4003: Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data;
4005: PetscInt i,m = A->rmap->n;
4008: c = (Mat_SeqAIJ*)C->data;
4010: C->factortype = A->factortype;
4011: c->row = 0;
4012: c->col = 0;
4013: c->icol = 0;
4014: c->reallocs = 0;
4016: C->assembled = PETSC_TRUE;
4018: PetscLayoutReference(A->rmap,&C->rmap);
4019: PetscLayoutReference(A->cmap,&C->cmap);
4021: PetscMalloc2(m,&c->imax,m,&c->ilen);
4022: PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));
4023: for (i=0; i<m; i++) {
4024: c->imax[i] = a->imax[i];
4025: c->ilen[i] = a->ilen[i];
4026: }
4028: /* allocate the matrix space */
4029: if (mallocmatspace) {
4030: PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);
4031: PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));
4033: c->singlemalloc = PETSC_TRUE;
4035: PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));
4036: if (m > 0) {
4037: PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));
4038: if (cpvalues == MAT_COPY_VALUES) {
4039: PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));
4040: } else {
4041: PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));
4042: }
4043: }
4044: }
4046: c->ignorezeroentries = a->ignorezeroentries;
4047: c->roworiented = a->roworiented;
4048: c->nonew = a->nonew;
4049: if (a->diag) {
4050: PetscMalloc1(m+1,&c->diag);
4051: PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));
4052: for (i=0; i<m; i++) {
4053: c->diag[i] = a->diag[i];
4054: }
4055: } else c->diag = 0;
4057: c->solve_work = 0;
4058: c->saved_values = 0;
4059: c->idiag = 0;
4060: c->ssor_work = 0;
4061: c->keepnonzeropattern = a->keepnonzeropattern;
4062: c->free_a = PETSC_TRUE;
4063: c->free_ij = PETSC_TRUE;
4065: c->rmax = a->rmax;
4066: c->nz = a->nz;
4067: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
4068: C->preallocated = PETSC_TRUE;
4070: c->compressedrow.use = a->compressedrow.use;
4071: c->compressedrow.nrows = a->compressedrow.nrows;
4072: if (a->compressedrow.use) {
4073: i = a->compressedrow.nrows;
4074: PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);
4075: PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));
4076: PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));
4077: } else {
4078: c->compressedrow.use = PETSC_FALSE;
4079: c->compressedrow.i = NULL;
4080: c->compressedrow.rindex = NULL;
4081: }
4082: c->nonzerorowcnt = a->nonzerorowcnt;
4083: C->nonzerostate = A->nonzerostate;
4085: MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);
4086: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
4087: return(0);
4088: }
4092: PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
4093: {
4097: MatCreate(PetscObjectComm((PetscObject)A),B);
4098: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
4099: if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) {
4100: MatSetBlockSizesFromMats(*B,A,A);
4101: }
4102: MatSetType(*B,((PetscObject)A)->type_name);
4103: MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);
4104: return(0);
4105: }
4109: PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer)
4110: {
4111: Mat_SeqAIJ *a;
4113: PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols;
4114: int fd;
4115: PetscMPIInt size;
4116: MPI_Comm comm;
4117: PetscInt bs = newMat->rmap->bs;
4120: /* force binary viewer to load .info file if it has not yet done so */
4121: PetscViewerSetUp(viewer);
4122: PetscObjectGetComm((PetscObject)viewer,&comm);
4123: MPI_Comm_size(comm,&size);
4124: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor");
4126: PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");
4127: PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);
4128: PetscOptionsEnd();
4129: if (bs < 0) bs = 1;
4130: MatSetBlockSize(newMat,bs);
4132: PetscViewerBinaryGetDescriptor(viewer,&fd);
4133: PetscBinaryRead(fd,header,4,PETSC_INT);
4134: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file");
4135: M = header[1]; N = header[2]; nz = header[3];
4137: if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ");
4139: /* read in row lengths */
4140: PetscMalloc1(M,&rowlengths);
4141: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
4143: /* check if sum of rowlengths is same as nz */
4144: for (i=0,sum=0; i< M; i++) sum +=rowlengths[i];
4145: if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %dD, sum-row-lengths = %D\n",nz,sum);
4147: /* set global size if not set already*/
4148: if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) {
4149: MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);
4150: } else {
4151: /* if sizes and type are already set, check if the matrix global sizes are correct */
4152: MatGetSize(newMat,&rows,&cols);
4153: if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */
4154: MatGetLocalSize(newMat,&rows,&cols);
4155: }
4156: 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);
4157: }
4158: MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);
4159: a = (Mat_SeqAIJ*)newMat->data;
4161: PetscBinaryRead(fd,a->j,nz,PETSC_INT);
4163: /* read in nonzero values */
4164: PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);
4166: /* set matrix "i" values */
4167: a->i[0] = 0;
4168: for (i=1; i<= M; i++) {
4169: a->i[i] = a->i[i-1] + rowlengths[i-1];
4170: a->ilen[i-1] = rowlengths[i-1];
4171: }
4172: PetscFree(rowlengths);
4174: MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);
4175: MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);
4176: return(0);
4177: }
4181: PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg)
4182: {
4183: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data;
4185: #if defined(PETSC_USE_COMPLEX)
4186: PetscInt k;
4187: #endif
4190: /* If the matrix dimensions are not equal,or no of nonzeros */
4191: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) {
4192: *flg = PETSC_FALSE;
4193: return(0);
4194: }
4196: /* if the a->i are the same */
4197: PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);
4198: if (!*flg) return(0);
4200: /* if a->j are the same */
4201: PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);
4202: if (!*flg) return(0);
4204: /* if a->a are the same */
4205: #if defined(PETSC_USE_COMPLEX)
4206: for (k=0; k<a->nz; k++) {
4207: if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) {
4208: *flg = PETSC_FALSE;
4209: return(0);
4210: }
4211: }
4212: #else
4213: PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);
4214: #endif
4215: return(0);
4216: }
4220: /*@
4221: MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format)
4222: provided by the user.
4224: Collective on MPI_Comm
4226: Input Parameters:
4227: + comm - must be an MPI communicator of size 1
4228: . m - number of rows
4229: . n - number of columns
4230: . i - row indices
4231: . j - column indices
4232: - a - matrix values
4234: Output Parameter:
4235: . mat - the matrix
4237: Level: intermediate
4239: Notes:
4240: The i, j, and a arrays are not copied by this routine, the user must free these arrays
4241: once the matrix is destroyed and not before
4243: You cannot set new nonzero locations into this matrix, that will generate an error.
4245: The i and j indices are 0 based
4247: The format which is used for the sparse matrix input, is equivalent to a
4248: row-major ordering.. i.e for the following matrix, the input data expected is
4249: as shown:
4251: 1 0 0
4252: 2 0 3
4253: 4 5 6
4255: i = {0,1,3,6} [size = nrow+1 = 3+1]
4256: j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row
4257: v = {1,2,3,4,5,6} [size = nz = 6]
4260: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
4262: @*/
4263: PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat)
4264: {
4266: PetscInt ii;
4267: Mat_SeqAIJ *aij;
4268: #if defined(PETSC_USE_DEBUG)
4269: PetscInt jj;
4270: #endif
4273: if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
4274: MatCreate(comm,mat);
4275: MatSetSizes(*mat,m,n,m,n);
4276: /* MatSetBlockSizes(*mat,,); */
4277: MatSetType(*mat,MATSEQAIJ);
4278: MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);
4279: aij = (Mat_SeqAIJ*)(*mat)->data;
4280: PetscMalloc2(m,&aij->imax,m,&aij->ilen);
4282: aij->i = i;
4283: aij->j = j;
4284: aij->a = a;
4285: aij->singlemalloc = PETSC_FALSE;
4286: aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
4287: aij->free_a = PETSC_FALSE;
4288: aij->free_ij = PETSC_FALSE;
4290: for (ii=0; ii<m; ii++) {
4291: aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
4292: #if defined(PETSC_USE_DEBUG)
4293: 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]);
4294: for (jj=i[ii]+1; jj<i[ii+1]; jj++) {
4295: if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii);
4296: if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii);
4297: }
4298: #endif
4299: }
4300: #if defined(PETSC_USE_DEBUG)
4301: for (ii=0; ii<aij->i[m]; ii++) {
4302: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]);
4303: 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]);
4304: }
4305: #endif
4307: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
4308: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
4309: return(0);
4310: }
4313: /*@C
4314: MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format)
4315: provided by the user.
4317: Collective on MPI_Comm
4319: Input Parameters:
4320: + comm - must be an MPI communicator of size 1
4321: . m - number of rows
4322: . n - number of columns
4323: . i - row indices
4324: . j - column indices
4325: . a - matrix values
4326: . nz - number of nonzeros
4327: - idx - 0 or 1 based
4329: Output Parameter:
4330: . mat - the matrix
4332: Level: intermediate
4334: Notes:
4335: The i and j indices are 0 based
4337: The format which is used for the sparse matrix input, is equivalent to a
4338: row-major ordering.. i.e for the following matrix, the input data expected is
4339: as shown:
4341: 1 0 0
4342: 2 0 3
4343: 4 5 6
4345: i = {0,1,1,2,2,2}
4346: j = {0,0,2,0,1,2}
4347: v = {1,2,3,4,5,6}
4350: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
4352: @*/
4353: PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx)
4354: {
4356: PetscInt ii, *nnz, one = 1,row,col;
4360: PetscCalloc1(m,&nnz);
4361: for (ii = 0; ii < nz; ii++) {
4362: nnz[i[ii] - !!idx] += 1;
4363: }
4364: MatCreate(comm,mat);
4365: MatSetSizes(*mat,m,n,m,n);
4366: MatSetType(*mat,MATSEQAIJ);
4367: MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);
4368: for (ii = 0; ii < nz; ii++) {
4369: if (idx) {
4370: row = i[ii] - 1;
4371: col = j[ii] - 1;
4372: } else {
4373: row = i[ii];
4374: col = j[ii];
4375: }
4376: MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);
4377: }
4378: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
4379: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
4380: PetscFree(nnz);
4381: return(0);
4382: }
4386: PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring)
4387: {
4389: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4392: if (coloring->ctype == IS_COLORING_GLOBAL) {
4393: ISColoringReference(coloring);
4394: a->coloring = coloring;
4395: } else if (coloring->ctype == IS_COLORING_GHOSTED) {
4396: PetscInt i,*larray;
4397: ISColoring ocoloring;
4398: ISColoringValue *colors;
4400: /* set coloring for diagonal portion */
4401: PetscMalloc1(A->cmap->n,&larray);
4402: for (i=0; i<A->cmap->n; i++) larray[i] = i;
4403: ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);
4404: PetscMalloc1(A->cmap->n,&colors);
4405: for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]];
4406: PetscFree(larray);
4407: ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,PETSC_OWN_POINTER,&ocoloring);
4408: a->coloring = ocoloring;
4409: }
4410: return(0);
4411: }
4415: PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues)
4416: {
4417: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4418: PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j;
4419: MatScalar *v = a->a;
4420: PetscScalar *values = (PetscScalar*)advalues;
4421: ISColoringValue *color;
4424: if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix");
4425: color = a->coloring->colors;
4426: /* loop over rows */
4427: for (i=0; i<m; i++) {
4428: nz = ii[i+1] - ii[i];
4429: /* loop over columns putting computed value into matrix */
4430: for (j=0; j<nz; j++) *v++ = values[color[*jj++]];
4431: values += nl; /* jump to next row of derivatives */
4432: }
4433: return(0);
4434: }
4438: PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A)
4439: {
4440: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
4444: a->idiagvalid = PETSC_FALSE;
4445: a->ibdiagvalid = PETSC_FALSE;
4447: MatSeqAIJInvalidateDiagonal_Inode(A);
4448: return(0);
4449: }
4453: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
4454: {
4458: MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);
4459: return(0);
4460: }
4462: /*
4463: Permute A into C's *local* index space using rowemb,colemb.
4464: The embedding are supposed to be injections and the above implies that the range of rowemb is a subset
4465: of [0,m), colemb is in [0,n).
4466: If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A.
4467: */
4470: PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B)
4471: {
4472: /* If making this function public, change the error returned in this function away from _PLIB. */
4474: Mat_SeqAIJ *Baij;
4475: PetscBool seqaij;
4476: PetscInt m,n,*nz,i,j,count;
4477: PetscScalar v;
4478: const PetscInt *rowindices,*colindices;
4481: if (!B) return(0);
4482: /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */
4483: PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);
4484: if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type");
4485: if (rowemb) {
4486: ISGetLocalSize(rowemb,&m);
4487: 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);
4488: } else {
4489: if (C->rmap->n != B->rmap->n) {
4490: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix");
4491: }
4492: }
4493: if (colemb) {
4494: ISGetLocalSize(colemb,&n);
4495: 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);
4496: } else {
4497: if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix");
4498: }
4500: Baij = (Mat_SeqAIJ*)(B->data);
4501: if (pattern == DIFFERENT_NONZERO_PATTERN) {
4502: PetscMalloc1(B->rmap->n,&nz);
4503: for (i=0; i<B->rmap->n; i++) {
4504: nz[i] = Baij->i[i+1] - Baij->i[i];
4505: }
4506: MatSeqAIJSetPreallocation(C,0,nz);
4507: PetscFree(nz);
4508: }
4509: if (pattern == SUBSET_NONZERO_PATTERN) {
4510: MatZeroEntries(C);
4511: }
4512: count = 0;
4513: rowindices = NULL;
4514: colindices = NULL;
4515: if (rowemb) {
4516: ISGetIndices(rowemb,&rowindices);
4517: }
4518: if (colemb) {
4519: ISGetIndices(colemb,&colindices);
4520: }
4521: for (i=0; i<B->rmap->n; i++) {
4522: PetscInt row;
4523: row = i;
4524: if (rowindices) row = rowindices[i];
4525: for (j=Baij->i[i]; j<Baij->i[i+1]; j++) {
4526: PetscInt col;
4527: col = Baij->j[count];
4528: if (colindices) col = colindices[col];
4529: v = Baij->a[count];
4530: MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);
4531: ++count;
4532: }
4533: }
4534: /* FIXME: set C's nonzerostate correctly. */
4535: /* Assembly for C is necessary. */
4536: C->preallocated = PETSC_TRUE;
4537: C->assembled = PETSC_TRUE;
4538: C->was_assembled = PETSC_FALSE;
4539: return(0);
4540: }
4543: /*
4544: Special version for direct calls from Fortran
4545: */
4546: #include <petsc/private/fortranimpl.h>
4547: #if defined(PETSC_HAVE_FORTRAN_CAPS)
4548: #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ
4549: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
4550: #define matsetvaluesseqaij_ matsetvaluesseqaij
4551: #endif
4553: /* Change these macros so can be used in void function */
4554: #undef CHKERRQ
4555: #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr)
4556: #undef SETERRQ2
4557: #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr)
4558: #undef SETERRQ3
4559: #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr)
4563: PETSC_EXTERN void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr)
4564: {
4565: Mat A = *AA;
4566: PetscInt m = *mm, n = *nn;
4567: InsertMode is = *isis;
4568: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4569: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
4570: PetscInt *imax,*ai,*ailen;
4572: PetscInt *aj,nonew = a->nonew,lastcol = -1;
4573: MatScalar *ap,value,*aa;
4574: PetscBool ignorezeroentries = a->ignorezeroentries;
4575: PetscBool roworiented = a->roworiented;
4578: MatCheckPreallocated(A,1);
4579: imax = a->imax;
4580: ai = a->i;
4581: ailen = a->ilen;
4582: aj = a->j;
4583: aa = a->a;
4585: for (k=0; k<m; k++) { /* loop over added rows */
4586: row = im[k];
4587: if (row < 0) continue;
4588: #if defined(PETSC_USE_DEBUG)
4589: if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
4590: #endif
4591: rp = aj + ai[row]; ap = aa + ai[row];
4592: rmax = imax[row]; nrow = ailen[row];
4593: low = 0;
4594: high = nrow;
4595: for (l=0; l<n; l++) { /* loop over added columns */
4596: if (in[l] < 0) continue;
4597: #if defined(PETSC_USE_DEBUG)
4598: if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
4599: #endif
4600: col = in[l];
4601: if (roworiented) value = v[l + k*n];
4602: else value = v[k + l*m];
4604: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
4606: if (col <= lastcol) low = 0;
4607: else high = nrow;
4608: lastcol = col;
4609: while (high-low > 5) {
4610: t = (low+high)/2;
4611: if (rp[t] > col) high = t;
4612: else low = t;
4613: }
4614: for (i=low; i<high; i++) {
4615: if (rp[i] > col) break;
4616: if (rp[i] == col) {
4617: if (is == ADD_VALUES) ap[i] += value;
4618: else ap[i] = value;
4619: goto noinsert;
4620: }
4621: }
4622: if (value == 0.0 && ignorezeroentries) goto noinsert;
4623: if (nonew == 1) goto noinsert;
4624: if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
4625: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
4626: N = nrow++ - 1; a->nz++; high++;
4627: /* shift up all the later entries in this row */
4628: for (ii=N; ii>=i; ii--) {
4629: rp[ii+1] = rp[ii];
4630: ap[ii+1] = ap[ii];
4631: }
4632: rp[i] = col;
4633: ap[i] = value;
4634: A->nonzerostate++;
4635: noinsert:;
4636: low = i + 1;
4637: }
4638: ailen[row] = nrow;
4639: }
4640: PetscFunctionReturnVoid();
4641: }