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