1
2 /*
3 Defines the basic matrix operations for the SBAIJ (compressed row)
4 matrix storage format.
5 */
6 #include <../src/mat/impls/baij/seq/baij.h> /*I "petscmat.h" I*/
7 #include <../src/mat/impls/sbaij/seq/sbaij.h>
8 #include <petscblaslapack.h>
9
10 #include <../src/mat/impls/sbaij/seq/relax.h>
11 #define USESHORT
12 #include <../src/mat/impls/sbaij/seq/relax.h>
13
14 #if defined(PETSC_HAVE_ELEMENTAL)
15 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat,MatType,MatReuse,Mat*);
16 #endif
17 #if defined(PETSC_HAVE_SCALAPACK)
18 PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*);
19 #endif
20 PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat,MatType,MatReuse,Mat*);
21
22 /*
23 Checks for missing diagonals
24 */
MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool * missing,PetscInt * dd)25 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd)
26 {
27 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
28 PetscErrorCode ierr;
29 PetscInt *diag,*ii = a->i,i;
30
31 PetscFunctionBegin;
32 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr);
33 *missing = PETSC_FALSE;
34 if (A->rmap->n > 0 && !ii) {
35 *missing = PETSC_TRUE;
36 if (dd) *dd = 0;
37 ierr = PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");CHKERRQ(ierr);
38 } else {
39 diag = a->diag;
40 for (i=0; i<a->mbs; i++) {
41 if (diag[i] >= ii[i+1]) {
42 *missing = PETSC_TRUE;
43 if (dd) *dd = i;
44 break;
45 }
46 }
47 }
48 PetscFunctionReturn(0);
49 }
50
MatMarkDiagonal_SeqSBAIJ(Mat A)51 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A)
52 {
53 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
54 PetscErrorCode ierr;
55 PetscInt i,j;
56
57 PetscFunctionBegin;
58 if (!a->diag) {
59 ierr = PetscMalloc1(a->mbs,&a->diag);CHKERRQ(ierr);
60 ierr = PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt));CHKERRQ(ierr);
61 a->free_diag = PETSC_TRUE;
62 }
63 for (i=0; i<a->mbs; i++) {
64 a->diag[i] = a->i[i+1];
65 for (j=a->i[i]; j<a->i[i+1]; j++) {
66 if (a->j[j] == i) {
67 a->diag[i] = j;
68 break;
69 }
70 }
71 }
72 PetscFunctionReturn(0);
73 }
74
MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt * nn,const PetscInt * inia[],const PetscInt * inja[],PetscBool * done)75 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done)
76 {
77 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
78 PetscErrorCode ierr;
79 PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt;
80 PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja;
81
82 PetscFunctionBegin;
83 *nn = n;
84 if (!ia) PetscFunctionReturn(0);
85 if (symmetric) {
86 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja);CHKERRQ(ierr);
87 nz = tia[n];
88 } else {
89 tia = a->i; tja = a->j;
90 }
91
92 if (!blockcompressed && bs > 1) {
93 (*nn) *= bs;
94 /* malloc & create the natural set of indices */
95 ierr = PetscMalloc1((n+1)*bs,ia);CHKERRQ(ierr);
96 if (n) {
97 (*ia)[0] = oshift;
98 for (j=1; j<bs; j++) {
99 (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1];
100 }
101 }
102
103 for (i=1; i<n; i++) {
104 (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1];
105 for (j=1; j<bs; j++) {
106 (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1];
107 }
108 }
109 if (n) {
110 (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1];
111 }
112
113 if (inja) {
114 ierr = PetscMalloc1(nz*bs*bs,ja);CHKERRQ(ierr);
115 cnt = 0;
116 for (i=0; i<n; i++) {
117 for (j=0; j<bs; j++) {
118 for (k=tia[i]; k<tia[i+1]; k++) {
119 for (l=0; l<bs; l++) {
120 (*ja)[cnt++] = bs*tja[k] + l;
121 }
122 }
123 }
124 }
125 }
126
127 if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
128 ierr = PetscFree(tia);CHKERRQ(ierr);
129 ierr = PetscFree(tja);CHKERRQ(ierr);
130 }
131 } else if (oshift == 1) {
132 if (symmetric) {
133 nz = tia[A->rmap->n/bs];
134 /* add 1 to i and j indices */
135 for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1;
136 *ia = tia;
137 if (ja) {
138 for (i=0; i<nz; i++) tja[i] = tja[i] + 1;
139 *ja = tja;
140 }
141 } else {
142 nz = a->i[A->rmap->n/bs];
143 /* malloc space and add 1 to i and j indices */
144 ierr = PetscMalloc1(A->rmap->n/bs+1,ia);CHKERRQ(ierr);
145 for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1;
146 if (ja) {
147 ierr = PetscMalloc1(nz,ja);CHKERRQ(ierr);
148 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1;
149 }
150 }
151 } else {
152 *ia = tia;
153 if (ja) *ja = tja;
154 }
155 PetscFunctionReturn(0);
156 }
157
MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt * nn,const PetscInt * ia[],const PetscInt * ja[],PetscBool * done)158 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
159 {
160 PetscErrorCode ierr;
161
162 PetscFunctionBegin;
163 if (!ia) PetscFunctionReturn(0);
164 if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
165 ierr = PetscFree(*ia);CHKERRQ(ierr);
166 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);}
167 }
168 PetscFunctionReturn(0);
169 }
170
MatDestroy_SeqSBAIJ(Mat A)171 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A)
172 {
173 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
174 PetscErrorCode ierr;
175
176 PetscFunctionBegin;
177 #if defined(PETSC_USE_LOG)
178 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz);
179 #endif
180 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr);
181 if (a->free_diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);}
182 ierr = ISDestroy(&a->row);CHKERRQ(ierr);
183 ierr = ISDestroy(&a->col);CHKERRQ(ierr);
184 ierr = ISDestroy(&a->icol);CHKERRQ(ierr);
185 ierr = PetscFree(a->idiag);CHKERRQ(ierr);
186 ierr = PetscFree(a->inode.size);CHKERRQ(ierr);
187 if (a->free_imax_ilen) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);}
188 ierr = PetscFree(a->solve_work);CHKERRQ(ierr);
189 ierr = PetscFree(a->sor_work);CHKERRQ(ierr);
190 ierr = PetscFree(a->solves_work);CHKERRQ(ierr);
191 ierr = PetscFree(a->mult_work);CHKERRQ(ierr);
192 ierr = PetscFree(a->saved_values);CHKERRQ(ierr);
193 if (a->free_jshort) {ierr = PetscFree(a->jshort);CHKERRQ(ierr);}
194 ierr = PetscFree(a->inew);CHKERRQ(ierr);
195 ierr = MatDestroy(&a->parent);CHKERRQ(ierr);
196 ierr = PetscFree(A->data);CHKERRQ(ierr);
197
198 ierr = PetscObjectChangeTypeName((PetscObject)A,NULL);CHKERRQ(ierr);
199 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr);
200 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr);
201 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL);CHKERRQ(ierr);
202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL);CHKERRQ(ierr);
203 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL);CHKERRQ(ierr);
204 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL);CHKERRQ(ierr);
205 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr);
206 #if defined(PETSC_HAVE_ELEMENTAL)
207 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL);CHKERRQ(ierr);
208 #endif
209 #if defined(PETSC_HAVE_SCALAPACK)
210 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_scalapack_C",NULL);CHKERRQ(ierr);
211 #endif
212 PetscFunctionReturn(0);
213 }
214
MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg)215 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg)
216 {
217 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
218 #if defined(PETSC_USE_COMPLEX)
219 PetscInt bs;
220 #endif
221 PetscErrorCode ierr;
222
223 PetscFunctionBegin;
224 #if defined(PETSC_USE_COMPLEX)
225 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
226 #endif
227 switch (op) {
228 case MAT_ROW_ORIENTED:
229 a->roworiented = flg;
230 break;
231 case MAT_KEEP_NONZERO_PATTERN:
232 a->keepnonzeropattern = flg;
233 break;
234 case MAT_NEW_NONZERO_LOCATIONS:
235 a->nonew = (flg ? 0 : 1);
236 break;
237 case MAT_NEW_NONZERO_LOCATION_ERR:
238 a->nonew = (flg ? -1 : 0);
239 break;
240 case MAT_NEW_NONZERO_ALLOCATION_ERR:
241 a->nonew = (flg ? -2 : 0);
242 break;
243 case MAT_UNUSED_NONZERO_LOCATION_ERR:
244 a->nounused = (flg ? -1 : 0);
245 break;
246 case MAT_NEW_DIAGONALS:
247 case MAT_IGNORE_OFF_PROC_ENTRIES:
248 case MAT_USE_HASH_TABLE:
249 case MAT_SORTED_FULL:
250 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr);
251 break;
252 case MAT_HERMITIAN:
253 #if defined(PETSC_USE_COMPLEX)
254 if (flg) { /* disable transpose ops */
255 if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1");
256 A->ops->multtranspose = NULL;
257 A->ops->multtransposeadd = NULL;
258 A->symmetric = PETSC_FALSE;
259 }
260 #endif
261 break;
262 case MAT_SYMMETRIC:
263 case MAT_SPD:
264 #if defined(PETSC_USE_COMPLEX)
265 if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */
266 A->ops->multtranspose = A->ops->mult;
267 A->ops->multtransposeadd = A->ops->multadd;
268 }
269 #endif
270 break;
271 /* These options are handled directly by MatSetOption() */
272 case MAT_STRUCTURALLY_SYMMETRIC:
273 case MAT_SYMMETRY_ETERNAL:
274 case MAT_STRUCTURE_ONLY:
275 /* These options are handled directly by MatSetOption() */
276 break;
277 case MAT_IGNORE_LOWER_TRIANGULAR:
278 a->ignore_ltriangular = flg;
279 break;
280 case MAT_ERROR_LOWER_TRIANGULAR:
281 a->ignore_ltriangular = flg;
282 break;
283 case MAT_GETROW_UPPERTRIANGULAR:
284 a->getrow_utriangular = flg;
285 break;
286 case MAT_SUBMAT_SINGLEIS:
287 break;
288 default:
289 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
290 }
291 PetscFunctionReturn(0);
292 }
293
MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt * nz,PetscInt ** idx,PetscScalar ** v)294 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
295 {
296 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
297 PetscErrorCode ierr;
298
299 PetscFunctionBegin;
300 if (A && !a->getrow_utriangular) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()");
301
302 /* Get the upper triangular part of the row */
303 ierr = MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);CHKERRQ(ierr);
304 PetscFunctionReturn(0);
305 }
306
MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt * nz,PetscInt ** idx,PetscScalar ** v)307 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
308 {
309 PetscErrorCode ierr;
310
311 PetscFunctionBegin;
312 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}
313 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);}
314 PetscFunctionReturn(0);
315 }
316
MatGetRowUpperTriangular_SeqSBAIJ(Mat A)317 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A)
318 {
319 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
320
321 PetscFunctionBegin;
322 a->getrow_utriangular = PETSC_TRUE;
323 PetscFunctionReturn(0);
324 }
325
MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)326 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A)
327 {
328 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
329
330 PetscFunctionBegin;
331 a->getrow_utriangular = PETSC_FALSE;
332 PetscFunctionReturn(0);
333 }
334
MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat * B)335 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B)
336 {
337 PetscErrorCode ierr;
338
339 PetscFunctionBegin;
340 if (reuse == MAT_INITIAL_MATRIX) {
341 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr);
342 } else if (reuse == MAT_REUSE_MATRIX) {
343 ierr = MatCopy(A,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
344 }
345 PetscFunctionReturn(0);
346 }
347
MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)348 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer)
349 {
350 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
351 PetscErrorCode ierr;
352 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
353 PetscViewerFormat format;
354 PetscInt *diag;
355
356 PetscFunctionBegin;
357 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
358 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
359 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr);
360 } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
361 Mat aij;
362 const char *matname;
363
364 if (A->factortype && bs>1) {
365 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");CHKERRQ(ierr);
366 PetscFunctionReturn(0);
367 }
368 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr);
369 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr);
370 ierr = PetscObjectSetName((PetscObject)aij,matname);CHKERRQ(ierr);
371 ierr = MatView(aij,viewer);CHKERRQ(ierr);
372 ierr = MatDestroy(&aij);CHKERRQ(ierr);
373 } else if (format == PETSC_VIEWER_ASCII_COMMON) {
374 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
375 for (i=0; i<a->mbs; i++) {
376 for (j=0; j<bs; j++) {
377 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr);
378 for (k=a->i[i]; k<a->i[i+1]; k++) {
379 for (l=0; l<bs; l++) {
380 #if defined(PETSC_USE_COMPLEX)
381 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
382 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
383 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
384 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
385 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
386 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
387 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
388 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
389 }
390 #else
391 if (a->a[bs2*k + l*bs + j] != 0.0) {
392 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr);
393 }
394 #endif
395 }
396 }
397 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
398 }
399 }
400 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
401 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
402 PetscFunctionReturn(0);
403 } else {
404 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
405 if (A->factortype) { /* for factored matrix */
406 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet");
407
408 diag=a->diag;
409 for (i=0; i<a->mbs; i++) { /* for row block i */
410 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr);
411 /* diagonal entry */
412 #if defined(PETSC_USE_COMPLEX)
413 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) {
414 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),(double)PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr);
415 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) {
416 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),-(double)PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr);
417 } else {
418 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr);
419 }
420 #else
421 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));CHKERRQ(ierr);
422 #endif
423 /* off-diagonal entries */
424 for (k=a->i[i]; k<a->i[i+1]-1; k++) {
425 #if defined(PETSC_USE_COMPLEX)
426 if (PetscImaginaryPart(a->a[k]) > 0.0) {
427 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr);
428 } else if (PetscImaginaryPart(a->a[k]) < 0.0) {
429 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr);
430 } else {
431 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));CHKERRQ(ierr);
432 }
433 #else
434 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);CHKERRQ(ierr);
435 #endif
436 }
437 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
438 }
439
440 } else { /* for non-factored matrix */
441 for (i=0; i<a->mbs; i++) { /* for row block i */
442 for (j=0; j<bs; j++) { /* for row bs*i + j */
443 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr);
444 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */
445 for (l=0; l<bs; l++) { /* for column */
446 #if defined(PETSC_USE_COMPLEX)
447 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
448 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
449 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
450 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
451 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
452 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
453 } else {
454 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr);
455 }
456 #else
457 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr);
458 #endif
459 }
460 }
461 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
462 }
463 }
464 }
465 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
466 }
467 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
468 PetscFunctionReturn(0);
469 }
470
471 #include <petscdraw.h>
MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void * Aa)472 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
473 {
474 Mat A = (Mat) Aa;
475 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data;
476 PetscErrorCode ierr;
477 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
478 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
479 MatScalar *aa;
480 PetscViewer viewer;
481
482 PetscFunctionBegin;
483 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr);
484 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr);
485
486 /* loop over matrix elements drawing boxes */
487
488 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr);
489 ierr = PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");CHKERRQ(ierr);
490 /* Blue for negative, Cyan for zero and Red for positive */
491 color = PETSC_DRAW_BLUE;
492 for (i=0,row=0; i<mbs; i++,row+=bs) {
493 for (j=a->i[i]; j<a->i[i+1]; j++) {
494 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
495 x_l = a->j[j]*bs; x_r = x_l + 1.0;
496 aa = a->a + j*bs2;
497 for (k=0; k<bs; k++) {
498 for (l=0; l<bs; l++) {
499 if (PetscRealPart(*aa++) >= 0.) continue;
500 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr);
501 }
502 }
503 }
504 }
505 color = PETSC_DRAW_CYAN;
506 for (i=0,row=0; i<mbs; i++,row+=bs) {
507 for (j=a->i[i]; j<a->i[i+1]; j++) {
508 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
509 x_l = a->j[j]*bs; x_r = x_l + 1.0;
510 aa = a->a + j*bs2;
511 for (k=0; k<bs; k++) {
512 for (l=0; l<bs; l++) {
513 if (PetscRealPart(*aa++) != 0.) continue;
514 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr);
515 }
516 }
517 }
518 }
519 color = PETSC_DRAW_RED;
520 for (i=0,row=0; i<mbs; i++,row+=bs) {
521 for (j=a->i[i]; j<a->i[i+1]; j++) {
522 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
523 x_l = a->j[j]*bs; x_r = x_l + 1.0;
524 aa = a->a + j*bs2;
525 for (k=0; k<bs; k++) {
526 for (l=0; l<bs; l++) {
527 if (PetscRealPart(*aa++) <= 0.) continue;
528 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr);
529 }
530 }
531 }
532 }
533 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr);
534 PetscFunctionReturn(0);
535 }
536
MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)537 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer)
538 {
539 PetscErrorCode ierr;
540 PetscReal xl,yl,xr,yr,w,h;
541 PetscDraw draw;
542 PetscBool isnull;
543
544 PetscFunctionBegin;
545 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
546 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr);
547 if (isnull) PetscFunctionReturn(0);
548
549 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
550 xr += w; yr += h; xl = -w; yl = -h;
551 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr);
552 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr);
553 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr);
554 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr);
555 ierr = PetscDrawSave(draw);CHKERRQ(ierr);
556 PetscFunctionReturn(0);
557 }
558
559 /* Used for both MPIBAIJ and MPISBAIJ matrices */
560 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary
561
MatView_SeqSBAIJ(Mat A,PetscViewer viewer)562 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer)
563 {
564 PetscErrorCode ierr;
565 PetscBool iascii,isbinary,isdraw;
566
567 PetscFunctionBegin;
568 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
569 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
571 if (iascii) {
572 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr);
573 } else if (isbinary) {
574 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr);
575 } else if (isdraw) {
576 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr);
577 } else {
578 Mat B;
579 const char *matname;
580 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
581 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr);
582 ierr = PetscObjectSetName((PetscObject)B,matname);CHKERRQ(ierr);
583 ierr = MatView(B,viewer);CHKERRQ(ierr);
584 ierr = MatDestroy(&B);CHKERRQ(ierr);
585 }
586 PetscFunctionReturn(0);
587 }
588
589
MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])590 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
591 {
592 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
593 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
594 PetscInt *ai = a->i,*ailen = a->ilen;
595 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
596 MatScalar *ap,*aa = a->a;
597
598 PetscFunctionBegin;
599 for (k=0; k<m; k++) { /* loop over rows */
600 row = im[k]; brow = row/bs;
601 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
602 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);
603 rp = aj + ai[brow]; ap = aa + bs2*ai[brow];
604 nrow = ailen[brow];
605 for (l=0; l<n; l++) { /* loop over columns */
606 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
607 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);
608 col = in[l];
609 bcol = col/bs;
610 cidx = col%bs;
611 ridx = row%bs;
612 high = nrow;
613 low = 0; /* assume unsorted */
614 while (high-low > 5) {
615 t = (low+high)/2;
616 if (rp[t] > bcol) high = t;
617 else low = t;
618 }
619 for (i=low; i<high; i++) {
620 if (rp[i] > bcol) break;
621 if (rp[i] == bcol) {
622 *v++ = ap[bs2*i+bs*cidx+ridx];
623 goto finished;
624 }
625 }
626 *v++ = 0.0;
627 finished:;
628 }
629 }
630 PetscFunctionReturn(0);
631 }
632
633
MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)634 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
635 {
636 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
637 PetscErrorCode ierr;
638 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
639 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen;
640 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
641 PetscBool roworiented=a->roworiented;
642 const PetscScalar *value = v;
643 MatScalar *ap,*aa = a->a,*bap;
644
645 PetscFunctionBegin;
646 if (roworiented) stepval = (n-1)*bs;
647 else stepval = (m-1)*bs;
648
649 for (k=0; k<m; k++) { /* loop over added rows */
650 row = im[k];
651 if (row < 0) continue;
652 if (PetscUnlikelyDebug(row >= a->mbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index row too large %D max %D",row,a->mbs-1);
653 rp = aj + ai[row];
654 ap = aa + bs2*ai[row];
655 rmax = imax[row];
656 nrow = ailen[row];
657 low = 0;
658 high = nrow;
659 for (l=0; l<n; l++) { /* loop over added columns */
660 if (in[l] < 0) continue;
661 col = in[l];
662 if (PetscUnlikelyDebug(col >= a->nbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index column too large %D max %D",col,a->nbs-1);
663 if (col < row) {
664 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */
665 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
666 }
667 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs;
668 else value = v + l*(stepval+bs)*bs + k*bs;
669
670 if (col <= lastcol) low = 0;
671 else high = nrow;
672
673 lastcol = col;
674 while (high-low > 7) {
675 t = (low+high)/2;
676 if (rp[t] > col) high = t;
677 else low = t;
678 }
679 for (i=low; i<high; i++) {
680 if (rp[i] > col) break;
681 if (rp[i] == col) {
682 bap = ap + bs2*i;
683 if (roworiented) {
684 if (is == ADD_VALUES) {
685 for (ii=0; ii<bs; ii++,value+=stepval) {
686 for (jj=ii; jj<bs2; jj+=bs) {
687 bap[jj] += *value++;
688 }
689 }
690 } else {
691 for (ii=0; ii<bs; ii++,value+=stepval) {
692 for (jj=ii; jj<bs2; jj+=bs) {
693 bap[jj] = *value++;
694 }
695 }
696 }
697 } else {
698 if (is == ADD_VALUES) {
699 for (ii=0; ii<bs; ii++,value+=stepval) {
700 for (jj=0; jj<bs; jj++) {
701 *bap++ += *value++;
702 }
703 }
704 } else {
705 for (ii=0; ii<bs; ii++,value+=stepval) {
706 for (jj=0; jj<bs; jj++) {
707 *bap++ = *value++;
708 }
709 }
710 }
711 }
712 goto noinsert2;
713 }
714 }
715 if (nonew == 1) goto noinsert2;
716 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new block index nonzero block (%D, %D) in the matrix", row, col);
717 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
718 N = nrow++ - 1; high++;
719 /* shift up all the later entries in this row */
720 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr);
721 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr);
722 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr);
723 rp[i] = col;
724 bap = ap + bs2*i;
725 if (roworiented) {
726 for (ii=0; ii<bs; ii++,value+=stepval) {
727 for (jj=ii; jj<bs2; jj+=bs) {
728 bap[jj] = *value++;
729 }
730 }
731 } else {
732 for (ii=0; ii<bs; ii++,value+=stepval) {
733 for (jj=0; jj<bs; jj++) {
734 *bap++ = *value++;
735 }
736 }
737 }
738 noinsert2:;
739 low = i;
740 }
741 ailen[row] = nrow;
742 }
743 PetscFunctionReturn(0);
744 }
745
746 /*
747 This is not yet used
748 */
MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)749 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A)
750 {
751 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
752 PetscErrorCode ierr;
753 const PetscInt *ai = a->i, *aj = a->j,*cols;
754 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts;
755 PetscBool flag;
756
757 PetscFunctionBegin;
758 ierr = PetscMalloc1(m,&ns);CHKERRQ(ierr);
759 while (i < m) {
760 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */
761 /* Limits the number of elements in a node to 'a->inode.limit' */
762 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
763 nzy = ai[j+1] - ai[j];
764 if (nzy != (nzx - j + i)) break;
765 ierr = PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);CHKERRQ(ierr);
766 if (!flag) break;
767 }
768 ns[node_count++] = blk_size;
769
770 i = j;
771 }
772 if (!a->inode.size && m && node_count > .9*m) {
773 ierr = PetscFree(ns);CHKERRQ(ierr);
774 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr);
775 } else {
776 a->inode.node_count = node_count;
777
778 ierr = PetscMalloc1(node_count,&a->inode.size);CHKERRQ(ierr);
779 ierr = PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));CHKERRQ(ierr);
780 ierr = PetscArraycpy(a->inode.size,ns,node_count);CHKERRQ(ierr);
781 ierr = PetscFree(ns);CHKERRQ(ierr);
782 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr);
783
784 /* count collections of adjacent columns in each inode */
785 row = 0;
786 cnt = 0;
787 for (i=0; i<node_count; i++) {
788 cols = aj + ai[row] + a->inode.size[i];
789 nz = ai[row+1] - ai[row] - a->inode.size[i];
790 for (j=1; j<nz; j++) {
791 if (cols[j] != cols[j-1]+1) cnt++;
792 }
793 cnt++;
794 row += a->inode.size[i];
795 }
796 ierr = PetscMalloc1(2*cnt,&counts);CHKERRQ(ierr);
797 cnt = 0;
798 row = 0;
799 for (i=0; i<node_count; i++) {
800 cols = aj + ai[row] + a->inode.size[i];
801 counts[2*cnt] = cols[0];
802 nz = ai[row+1] - ai[row] - a->inode.size[i];
803 cnt2 = 1;
804 for (j=1; j<nz; j++) {
805 if (cols[j] != cols[j-1]+1) {
806 counts[2*(cnt++)+1] = cnt2;
807 counts[2*cnt] = cols[j];
808 cnt2 = 1;
809 } else cnt2++;
810 }
811 counts[2*(cnt++)+1] = cnt2;
812 row += a->inode.size[i];
813 }
814 ierr = PetscIntView(2*cnt,counts,NULL);CHKERRQ(ierr);
815 }
816 PetscFunctionReturn(0);
817 }
818
MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)819 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode)
820 {
821 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
822 PetscErrorCode ierr;
823 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
824 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen;
825 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0;
826 MatScalar *aa = a->a,*ap;
827
828 PetscFunctionBegin;
829 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0);
830
831 if (m) rmax = ailen[0];
832 for (i=1; i<mbs; i++) {
833 /* move each row back by the amount of empty slots (fshift) before it*/
834 fshift += imax[i-1] - ailen[i-1];
835 rmax = PetscMax(rmax,ailen[i]);
836 if (fshift) {
837 ip = aj + ai[i];
838 ap = aa + bs2*ai[i];
839 N = ailen[i];
840 ierr = PetscArraymove(ip-fshift,ip,N);CHKERRQ(ierr);
841 ierr = PetscArraymove(ap-bs2*fshift,ap,bs2*N);CHKERRQ(ierr);
842 }
843 ai[i] = ai[i-1] + ailen[i-1];
844 }
845 if (mbs) {
846 fshift += imax[mbs-1] - ailen[mbs-1];
847 ai[mbs] = ai[mbs-1] + ailen[mbs-1];
848 }
849 /* reset ilen and imax for each row */
850 for (i=0; i<mbs; i++) {
851 ailen[i] = imax[i] = ai[i+1] - ai[i];
852 }
853 a->nz = ai[mbs];
854
855 /* diagonals may have moved, reset it */
856 if (a->diag) {
857 ierr = PetscArraycpy(a->diag,ai,mbs);CHKERRQ(ierr);
858 }
859 if (fshift && a->nounused == -1) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2);
860
861 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr);
862 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr);
863 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr);
864
865 A->info.mallocs += a->reallocs;
866 a->reallocs = 0;
867 A->info.nz_unneeded = (PetscReal)fshift*bs2;
868 a->idiagvalid = PETSC_FALSE;
869 a->rmax = rmax;
870
871 if (A->cmap->n < 65536 && A->cmap->bs == 1) {
872 if (a->jshort && a->free_jshort) {
873 /* when matrix data structure is changed, previous jshort must be replaced */
874 ierr = PetscFree(a->jshort);CHKERRQ(ierr);
875 }
876 ierr = PetscMalloc1(a->i[A->rmap->n],&a->jshort);CHKERRQ(ierr);
877 ierr = PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr);
878 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i];
879 A->ops->mult = MatMult_SeqSBAIJ_1_ushort;
880 A->ops->sor = MatSOR_SeqSBAIJ_ushort;
881 a->free_jshort = PETSC_TRUE;
882 }
883 PetscFunctionReturn(0);
884 }
885
886 /*
887 This function returns an array of flags which indicate the locations of contiguous
888 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
889 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
890 Assume: sizes should be long enough to hold all the values.
891 */
MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[],PetscInt * bs_max)892 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
893 {
894 PetscInt i,j,k,row;
895 PetscBool flg;
896
897 PetscFunctionBegin;
898 for (i=0,j=0; i<n; j++) {
899 row = idx[i];
900 if (row%bs!=0) { /* Not the begining of a block */
901 sizes[j] = 1;
902 i++;
903 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */
904 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */
905 i++;
906 } else { /* Begining of the block, so check if the complete block exists */
907 flg = PETSC_TRUE;
908 for (k=1; k<bs; k++) {
909 if (row+k != idx[i+k]) { /* break in the block */
910 flg = PETSC_FALSE;
911 break;
912 }
913 }
914 if (flg) { /* No break in the bs */
915 sizes[j] = bs;
916 i += bs;
917 } else {
918 sizes[j] = 1;
919 i++;
920 }
921 }
922 }
923 *bs_max = j;
924 PetscFunctionReturn(0);
925 }
926
927
928 /* Only add/insert a(i,j) with i<=j (blocks).
929 Any a(i,j) with i>j input by user is ingored.
930 */
931
MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)932 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
933 {
934 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
935 PetscErrorCode ierr;
936 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
937 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented;
938 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
939 PetscInt ridx,cidx,bs2=a->bs2;
940 MatScalar *ap,value,*aa=a->a,*bap;
941
942 PetscFunctionBegin;
943 for (k=0; k<m; k++) { /* loop over added rows */
944 row = im[k]; /* row number */
945 brow = row/bs; /* block row number */
946 if (row < 0) continue;
947 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);
948 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/
949 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/
950 rmax = imax[brow]; /* maximum space allocated for this row */
951 nrow = ailen[brow]; /* actual length of this row */
952 low = 0;
953 high = nrow;
954 for (l=0; l<n; l++) { /* loop over added columns */
955 if (in[l] < 0) continue;
956 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);
957 col = in[l];
958 bcol = col/bs; /* block col number */
959
960 if (brow > bcol) {
961 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */
962 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
963 }
964
965 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */
966 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) {
967 /* element value a(k,l) */
968 if (roworiented) value = v[l + k*n];
969 else value = v[k + l*m];
970
971 /* move pointer bap to a(k,l) quickly and add/insert value */
972 if (col <= lastcol) low = 0;
973 else high = nrow;
974
975 lastcol = col;
976 while (high-low > 7) {
977 t = (low+high)/2;
978 if (rp[t] > bcol) high = t;
979 else low = t;
980 }
981 for (i=low; i<high; i++) {
982 if (rp[i] > bcol) break;
983 if (rp[i] == bcol) {
984 bap = ap + bs2*i + bs*cidx + ridx;
985 if (is == ADD_VALUES) *bap += value;
986 else *bap = value;
987 /* for diag block, add/insert its symmetric element a(cidx,ridx) */
988 if (brow == bcol && ridx < cidx) {
989 bap = ap + bs2*i + bs*ridx + cidx;
990 if (is == ADD_VALUES) *bap += value;
991 else *bap = value;
992 }
993 goto noinsert1;
994 }
995 }
996
997 if (nonew == 1) goto noinsert1;
998 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
999 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
1000
1001 N = nrow++ - 1; high++;
1002 /* shift up all the later entries in this row */
1003 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr);
1004 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr);
1005 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr);
1006 rp[i] = bcol;
1007 ap[bs2*i + bs*cidx + ridx] = value;
1008 /* for diag block, add/insert its symmetric element a(cidx,ridx) */
1009 if (brow == bcol && ridx < cidx) {
1010 ap[bs2*i + bs*ridx + cidx] = value;
1011 }
1012 A->nonzerostate++;
1013 noinsert1:;
1014 low = i;
1015 }
1016 } /* end of loop over added columns */
1017 ailen[brow] = nrow;
1018 } /* end of loop over added rows */
1019 PetscFunctionReturn(0);
1020 }
1021
MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo * info)1022 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info)
1023 {
1024 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data;
1025 Mat outA;
1026 PetscErrorCode ierr;
1027 PetscBool row_identity;
1028
1029 PetscFunctionBegin;
1030 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc");
1031 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr);
1032 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported");
1033 if (inA->rmap->bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */
1034
1035 outA = inA;
1036 inA->factortype = MAT_FACTOR_ICC;
1037 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr);
1038 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr);
1039
1040 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr);
1041 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr);
1042
1043 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr);
1044 ierr = ISDestroy(&a->row);CHKERRQ(ierr);
1045 a->row = row;
1046 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr);
1047 ierr = ISDestroy(&a->col);CHKERRQ(ierr);
1048 a->col = row;
1049
1050 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */
1051 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);}
1052 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr);
1053
1054 if (!a->solve_work) {
1055 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr);
1056 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr);
1057 }
1058
1059 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr);
1060 PetscFunctionReturn(0);
1061 }
1062
MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt * indices)1063 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices)
1064 {
1065 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data;
1066 PetscInt i,nz,n;
1067 PetscErrorCode ierr;
1068
1069 PetscFunctionBegin;
1070 nz = baij->maxnz;
1071 n = mat->cmap->n;
1072 for (i=0; i<nz; i++) baij->j[i] = indices[i];
1073
1074 baij->nz = nz;
1075 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];
1076
1077 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1078 PetscFunctionReturn(0);
1079 }
1080
1081 /*@
1082 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows
1083 in the matrix.
1084
1085 Input Parameters:
1086 + mat - the SeqSBAIJ matrix
1087 - indices - the column indices
1088
1089 Level: advanced
1090
1091 Notes:
1092 This can be called if you have precomputed the nonzero structure of the
1093 matrix and want to provide it to the matrix object to improve the performance
1094 of the MatSetValues() operation.
1095
1096 You MUST have set the correct numbers of nonzeros per row in the call to
1097 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted.
1098
1099 MUST be called before any calls to MatSetValues()
1100
1101 .seealso: MatCreateSeqSBAIJ
1102 @*/
MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt * indices)1103 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices)
1104 {
1105 PetscErrorCode ierr;
1106
1107 PetscFunctionBegin;
1108 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
1109 PetscValidPointer(indices,2);
1110 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr);
1111 PetscFunctionReturn(0);
1112 }
1113
MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)1114 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str)
1115 {
1116 PetscErrorCode ierr;
1117 PetscBool isbaij;
1118
1119 PetscFunctionBegin;
1120 ierr = PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr);
1121 if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name);
1122 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */
1123 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1124 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1125 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1126
1127 if (a->i[a->mbs] != b->i[b->mbs]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different");
1128 if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different");
1129 if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size");
1130 ierr = PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);CHKERRQ(ierr);
1131 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr);
1132 } else {
1133 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr);
1134 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr);
1135 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr);
1136 }
1137 PetscFunctionReturn(0);
1138 }
1139
MatSetUp_SeqSBAIJ(Mat A)1140 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A)
1141 {
1142 PetscErrorCode ierr;
1143
1144 PetscFunctionBegin;
1145 ierr = MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);CHKERRQ(ierr);
1146 PetscFunctionReturn(0);
1147 }
1148
MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar * array[])1149 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1150 {
1151 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1152
1153 PetscFunctionBegin;
1154 *array = a->a;
1155 PetscFunctionReturn(0);
1156 }
1157
MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar * array[])1158 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[])
1159 {
1160 PetscFunctionBegin;
1161 *array = NULL;
1162 PetscFunctionReturn(0);
1163 }
1164
MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt * nnz)1165 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz)
1166 {
1167 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
1168 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data;
1169 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data;
1170 PetscErrorCode ierr;
1171
1172 PetscFunctionBegin;
1173 /* Set the number of nonzeros in the new matrix */
1174 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr);
1175 PetscFunctionReturn(0);
1176 }
1177
MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)1178 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
1179 {
1180 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data;
1181 PetscErrorCode ierr;
1182 PetscInt bs=Y->rmap->bs,bs2=bs*bs;
1183 PetscBLASInt one = 1;
1184
1185 PetscFunctionBegin;
1186 if (str == SAME_NONZERO_PATTERN) {
1187 PetscScalar alpha = a;
1188 PetscBLASInt bnz;
1189 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr);
1190 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
1191 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);
1192 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1193 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);
1194 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr);
1195 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr);
1196 } else {
1197 Mat B;
1198 PetscInt *nnz;
1199 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
1200 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr);
1201 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr);
1202 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr);
1203 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr);
1204 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr);
1205 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr);
1206 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr);
1207 ierr = MatSetType(B,((PetscObject)Y)->type_name);CHKERRQ(ierr);
1208 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr);
1209 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr);
1210
1211 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr);
1212
1213 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr);
1214 ierr = PetscFree(nnz);CHKERRQ(ierr);
1215 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr);
1216 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr);
1217 }
1218 PetscFunctionReturn(0);
1219 }
1220
MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool * flg)1221 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1222 {
1223 PetscFunctionBegin;
1224 *flg = PETSC_TRUE;
1225 PetscFunctionReturn(0);
1226 }
1227
MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool * flg)1228 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg)
1229 {
1230 PetscFunctionBegin;
1231 *flg = PETSC_TRUE;
1232 PetscFunctionReturn(0);
1233 }
1234
MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool * flg)1235 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg)
1236 {
1237 PetscFunctionBegin;
1238 *flg = PETSC_FALSE;
1239 PetscFunctionReturn(0);
1240 }
1241
MatRealPart_SeqSBAIJ(Mat A)1242 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A)
1243 {
1244 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1245 PetscInt i,nz = a->bs2*a->i[a->mbs];
1246 MatScalar *aa = a->a;
1247
1248 PetscFunctionBegin;
1249 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1250 PetscFunctionReturn(0);
1251 }
1252
MatImaginaryPart_SeqSBAIJ(Mat A)1253 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A)
1254 {
1255 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
1256 PetscInt i,nz = a->bs2*a->i[a->mbs];
1257 MatScalar *aa = a->a;
1258
1259 PetscFunctionBegin;
1260 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1261 PetscFunctionReturn(0);
1262 }
1263
MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x,Vec b)1264 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
1265 {
1266 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data;
1267 PetscErrorCode ierr;
1268 PetscInt i,j,k,count;
1269 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col;
1270 PetscScalar zero = 0.0;
1271 MatScalar *aa;
1272 const PetscScalar *xx;
1273 PetscScalar *bb;
1274 PetscBool *zeroed,vecs = PETSC_FALSE;
1275
1276 PetscFunctionBegin;
1277 /* fix right hand side if needed */
1278 if (x && b) {
1279 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
1280 ierr = VecGetArray(b,&bb);CHKERRQ(ierr);
1281 vecs = PETSC_TRUE;
1282 }
1283
1284 /* zero the columns */
1285 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr);
1286 for (i=0; i<is_n; i++) {
1287 if (is_idx[i] < 0 || is_idx[i] >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",is_idx[i]);
1288 zeroed[is_idx[i]] = PETSC_TRUE;
1289 }
1290 if (vecs) {
1291 for (i=0; i<A->rmap->N; i++) {
1292 row = i/bs;
1293 for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1294 for (k=0; k<bs; k++) {
1295 col = bs*baij->j[j] + k;
1296 if (col <= i) continue;
1297 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1298 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col];
1299 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i];
1300 }
1301 }
1302 }
1303 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]];
1304 }
1305
1306 for (i=0; i<A->rmap->N; i++) {
1307 if (!zeroed[i]) {
1308 row = i/bs;
1309 for (j=baij->i[row]; j<baij->i[row+1]; j++) {
1310 for (k=0; k<bs; k++) {
1311 col = bs*baij->j[j] + k;
1312 if (zeroed[col]) {
1313 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
1314 aa[0] = 0.0;
1315 }
1316 }
1317 }
1318 }
1319 }
1320 ierr = PetscFree(zeroed);CHKERRQ(ierr);
1321 if (vecs) {
1322 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
1323 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr);
1324 }
1325
1326 /* zero the rows */
1327 for (i=0; i<is_n; i++) {
1328 row = is_idx[i];
1329 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
1330 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
1331 for (k=0; k<count; k++) {
1332 aa[0] = zero;
1333 aa += bs;
1334 }
1335 if (diag != 0.0) {
1336 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr);
1337 }
1338 }
1339 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1340 PetscFunctionReturn(0);
1341 }
1342
MatShift_SeqSBAIJ(Mat Y,PetscScalar a)1343 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a)
1344 {
1345 PetscErrorCode ierr;
1346 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data;
1347
1348 PetscFunctionBegin;
1349 if (!Y->preallocated || !aij->nz) {
1350 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr);
1351 }
1352 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr);
1353 PetscFunctionReturn(0);
1354 }
1355
1356 /* -------------------------------------------------------------------*/
1357 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ,
1358 MatGetRow_SeqSBAIJ,
1359 MatRestoreRow_SeqSBAIJ,
1360 MatMult_SeqSBAIJ_N,
1361 /* 4*/ MatMultAdd_SeqSBAIJ_N,
1362 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */
1363 MatMultAdd_SeqSBAIJ_N,
1364 NULL,
1365 NULL,
1366 NULL,
1367 /* 10*/ NULL,
1368 NULL,
1369 MatCholeskyFactor_SeqSBAIJ,
1370 MatSOR_SeqSBAIJ,
1371 MatTranspose_SeqSBAIJ,
1372 /* 15*/ MatGetInfo_SeqSBAIJ,
1373 MatEqual_SeqSBAIJ,
1374 MatGetDiagonal_SeqSBAIJ,
1375 MatDiagonalScale_SeqSBAIJ,
1376 MatNorm_SeqSBAIJ,
1377 /* 20*/ NULL,
1378 MatAssemblyEnd_SeqSBAIJ,
1379 MatSetOption_SeqSBAIJ,
1380 MatZeroEntries_SeqSBAIJ,
1381 /* 24*/ NULL,
1382 NULL,
1383 NULL,
1384 NULL,
1385 NULL,
1386 /* 29*/ MatSetUp_SeqSBAIJ,
1387 NULL,
1388 NULL,
1389 NULL,
1390 NULL,
1391 /* 34*/ MatDuplicate_SeqSBAIJ,
1392 NULL,
1393 NULL,
1394 NULL,
1395 MatICCFactor_SeqSBAIJ,
1396 /* 39*/ MatAXPY_SeqSBAIJ,
1397 MatCreateSubMatrices_SeqSBAIJ,
1398 MatIncreaseOverlap_SeqSBAIJ,
1399 MatGetValues_SeqSBAIJ,
1400 MatCopy_SeqSBAIJ,
1401 /* 44*/ NULL,
1402 MatScale_SeqSBAIJ,
1403 MatShift_SeqSBAIJ,
1404 NULL,
1405 MatZeroRowsColumns_SeqSBAIJ,
1406 /* 49*/ NULL,
1407 MatGetRowIJ_SeqSBAIJ,
1408 MatRestoreRowIJ_SeqSBAIJ,
1409 NULL,
1410 NULL,
1411 /* 54*/ NULL,
1412 NULL,
1413 NULL,
1414 NULL,
1415 MatSetValuesBlocked_SeqSBAIJ,
1416 /* 59*/ MatCreateSubMatrix_SeqSBAIJ,
1417 NULL,
1418 NULL,
1419 NULL,
1420 NULL,
1421 /* 64*/ NULL,
1422 NULL,
1423 NULL,
1424 NULL,
1425 NULL,
1426 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ,
1427 NULL,
1428 MatConvert_MPISBAIJ_Basic,
1429 NULL,
1430 NULL,
1431 /* 74*/ NULL,
1432 NULL,
1433 NULL,
1434 NULL,
1435 NULL,
1436 /* 79*/ NULL,
1437 NULL,
1438 NULL,
1439 MatGetInertia_SeqSBAIJ,
1440 MatLoad_SeqSBAIJ,
1441 /* 84*/ MatIsSymmetric_SeqSBAIJ,
1442 MatIsHermitian_SeqSBAIJ,
1443 MatIsStructurallySymmetric_SeqSBAIJ,
1444 NULL,
1445 NULL,
1446 /* 89*/ NULL,
1447 NULL,
1448 NULL,
1449 NULL,
1450 NULL,
1451 /* 94*/ NULL,
1452 NULL,
1453 NULL,
1454 NULL,
1455 NULL,
1456 /* 99*/ NULL,
1457 NULL,
1458 NULL,
1459 NULL,
1460 NULL,
1461 /*104*/ NULL,
1462 MatRealPart_SeqSBAIJ,
1463 MatImaginaryPart_SeqSBAIJ,
1464 MatGetRowUpperTriangular_SeqSBAIJ,
1465 MatRestoreRowUpperTriangular_SeqSBAIJ,
1466 /*109*/ NULL,
1467 NULL,
1468 NULL,
1469 NULL,
1470 MatMissingDiagonal_SeqSBAIJ,
1471 /*114*/ NULL,
1472 NULL,
1473 NULL,
1474 NULL,
1475 NULL,
1476 /*119*/ NULL,
1477 NULL,
1478 NULL,
1479 NULL,
1480 NULL,
1481 /*124*/ NULL,
1482 NULL,
1483 NULL,
1484 NULL,
1485 NULL,
1486 /*129*/ NULL,
1487 NULL,
1488 NULL,
1489 NULL,
1490 NULL,
1491 /*134*/ NULL,
1492 NULL,
1493 NULL,
1494 NULL,
1495 NULL,
1496 /*139*/ MatSetBlockSizes_Default,
1497 NULL,
1498 NULL,
1499 NULL,
1500 NULL,
1501 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ
1502 };
1503
MatStoreValues_SeqSBAIJ(Mat mat)1504 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat)
1505 {
1506 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1507 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1508 PetscErrorCode ierr;
1509
1510 PetscFunctionBegin;
1511 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1512
1513 /* allocate space for values if not already there */
1514 if (!aij->saved_values) {
1515 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr);
1516 }
1517
1518 /* copy values over */
1519 ierr = PetscArraycpy(aij->saved_values,aij->a,nz);CHKERRQ(ierr);
1520 PetscFunctionReturn(0);
1521 }
1522
MatRetrieveValues_SeqSBAIJ(Mat mat)1523 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat)
1524 {
1525 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data;
1526 PetscErrorCode ierr;
1527 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2;
1528
1529 PetscFunctionBegin;
1530 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1531 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
1532
1533 /* copy values over */
1534 ierr = PetscArraycpy(aij->a,aij->saved_values,nz);CHKERRQ(ierr);
1535 PetscFunctionReturn(0);
1536 }
1537
MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt * nnz)1538 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
1539 {
1540 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data;
1541 PetscErrorCode ierr;
1542 PetscInt i,mbs,nbs,bs2;
1543 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE;
1544
1545 PetscFunctionBegin;
1546 if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
1547
1548 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr);
1549 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
1550 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
1551 if (B->rmap->N > B->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"SEQSBAIJ matrix cannot have more rows %D than columns %D",B->rmap->N,B->cmap->N);
1552 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr);
1553
1554 B->preallocated = PETSC_TRUE;
1555
1556 mbs = B->rmap->N/bs;
1557 nbs = B->cmap->n/bs;
1558 bs2 = bs*bs;
1559
1560 if (mbs*bs != B->rmap->N || nbs*bs!=B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize");
1561
1562 if (nz == MAT_SKIP_ALLOCATION) {
1563 skipallocation = PETSC_TRUE;
1564 nz = 0;
1565 }
1566
1567 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3;
1568 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
1569 if (nnz) {
1570 for (i=0; i<mbs; i++) {
1571 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]);
1572 if (nnz[i] > nbs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D block rowlength %D",i,nnz[i],nbs);
1573 }
1574 }
1575
1576 B->ops->mult = MatMult_SeqSBAIJ_N;
1577 B->ops->multadd = MatMultAdd_SeqSBAIJ_N;
1578 B->ops->multtranspose = MatMult_SeqSBAIJ_N;
1579 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N;
1580
1581 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr);
1582 if (!flg) {
1583 switch (bs) {
1584 case 1:
1585 B->ops->mult = MatMult_SeqSBAIJ_1;
1586 B->ops->multadd = MatMultAdd_SeqSBAIJ_1;
1587 B->ops->multtranspose = MatMult_SeqSBAIJ_1;
1588 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1;
1589 break;
1590 case 2:
1591 B->ops->mult = MatMult_SeqSBAIJ_2;
1592 B->ops->multadd = MatMultAdd_SeqSBAIJ_2;
1593 B->ops->multtranspose = MatMult_SeqSBAIJ_2;
1594 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2;
1595 break;
1596 case 3:
1597 B->ops->mult = MatMult_SeqSBAIJ_3;
1598 B->ops->multadd = MatMultAdd_SeqSBAIJ_3;
1599 B->ops->multtranspose = MatMult_SeqSBAIJ_3;
1600 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3;
1601 break;
1602 case 4:
1603 B->ops->mult = MatMult_SeqSBAIJ_4;
1604 B->ops->multadd = MatMultAdd_SeqSBAIJ_4;
1605 B->ops->multtranspose = MatMult_SeqSBAIJ_4;
1606 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4;
1607 break;
1608 case 5:
1609 B->ops->mult = MatMult_SeqSBAIJ_5;
1610 B->ops->multadd = MatMultAdd_SeqSBAIJ_5;
1611 B->ops->multtranspose = MatMult_SeqSBAIJ_5;
1612 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5;
1613 break;
1614 case 6:
1615 B->ops->mult = MatMult_SeqSBAIJ_6;
1616 B->ops->multadd = MatMultAdd_SeqSBAIJ_6;
1617 B->ops->multtranspose = MatMult_SeqSBAIJ_6;
1618 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6;
1619 break;
1620 case 7:
1621 B->ops->mult = MatMult_SeqSBAIJ_7;
1622 B->ops->multadd = MatMultAdd_SeqSBAIJ_7;
1623 B->ops->multtranspose = MatMult_SeqSBAIJ_7;
1624 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7;
1625 break;
1626 }
1627 }
1628
1629 b->mbs = mbs;
1630 b->nbs = nbs;
1631 if (!skipallocation) {
1632 if (!b->imax) {
1633 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr);
1634
1635 b->free_imax_ilen = PETSC_TRUE;
1636
1637 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr);
1638 }
1639 if (!nnz) {
1640 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
1641 else if (nz <= 0) nz = 1;
1642 nz = PetscMin(nbs,nz);
1643 for (i=0; i<mbs; i++) b->imax[i] = nz;
1644 nz = nz*mbs; /* total nz */
1645 } else {
1646 PetscInt64 nz64 = 0;
1647 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
1648 ierr = PetscIntCast(nz64,&nz);CHKERRQ(ierr);
1649 }
1650 /* b->ilen will count nonzeros in each block row so far. */
1651 for (i=0; i<mbs; i++) b->ilen[i] = 0;
1652 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */
1653
1654 /* allocate the matrix space */
1655 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr);
1656 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr);
1657 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr);
1658 ierr = PetscArrayzero(b->a,nz*bs2);CHKERRQ(ierr);
1659 ierr = PetscArrayzero(b->j,nz);CHKERRQ(ierr);
1660
1661 b->singlemalloc = PETSC_TRUE;
1662
1663 /* pointer to beginning of each row */
1664 b->i[0] = 0;
1665 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1];
1666
1667 b->free_a = PETSC_TRUE;
1668 b->free_ij = PETSC_TRUE;
1669 } else {
1670 b->free_a = PETSC_FALSE;
1671 b->free_ij = PETSC_FALSE;
1672 }
1673
1674 b->bs2 = bs2;
1675 b->nz = 0;
1676 b->maxnz = nz;
1677 b->inew = NULL;
1678 b->jnew = NULL;
1679 b->anew = NULL;
1680 b->a2anew = NULL;
1681 b->permute = PETSC_FALSE;
1682
1683 B->was_assembled = PETSC_FALSE;
1684 B->assembled = PETSC_FALSE;
1685 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);}
1686 PetscFunctionReturn(0);
1687 }
1688
MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[],const PetscScalar V[])1689 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[])
1690 {
1691 PetscInt i,j,m,nz,anz, nz_max=0,*nnz;
1692 PetscScalar *values=NULL;
1693 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented;
1694 PetscErrorCode ierr;
1695
1696 PetscFunctionBegin;
1697 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
1698 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr);
1699 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr);
1700 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
1701 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
1702 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr);
1703 m = B->rmap->n/bs;
1704
1705 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]);
1706 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr);
1707 for (i=0; i<m; i++) {
1708 nz = ii[i+1] - ii[i];
1709 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz);
1710 anz = 0;
1711 for (j=0; j<nz; j++) {
1712 /* count only values on the diagonal or above */
1713 if (jj[ii[i] + j] >= i) {
1714 anz = nz - j;
1715 break;
1716 }
1717 }
1718 nz_max = PetscMax(nz_max,anz);
1719 nnz[i] = anz;
1720 }
1721 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr);
1722 ierr = PetscFree(nnz);CHKERRQ(ierr);
1723
1724 values = (PetscScalar*)V;
1725 if (!values) {
1726 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr);
1727 }
1728 for (i=0; i<m; i++) {
1729 PetscInt ncols = ii[i+1] - ii[i];
1730 const PetscInt *icols = jj + ii[i];
1731 if (!roworiented || bs == 1) {
1732 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
1733 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr);
1734 } else {
1735 for (j=0; j<ncols; j++) {
1736 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
1737 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr);
1738 }
1739 }
1740 }
1741 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); }
1742 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1743 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1744 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1745 PetscFunctionReturn(0);
1746 }
1747
1748 /*
1749 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization
1750 */
MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)1751 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural)
1752 {
1753 PetscErrorCode ierr;
1754 PetscBool flg = PETSC_FALSE;
1755 PetscInt bs = B->rmap->bs;
1756
1757 PetscFunctionBegin;
1758 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr);
1759 if (flg) bs = 8;
1760
1761 if (!natural) {
1762 switch (bs) {
1763 case 1:
1764 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace;
1765 break;
1766 case 2:
1767 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2;
1768 break;
1769 case 3:
1770 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3;
1771 break;
1772 case 4:
1773 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4;
1774 break;
1775 case 5:
1776 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5;
1777 break;
1778 case 6:
1779 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6;
1780 break;
1781 case 7:
1782 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7;
1783 break;
1784 default:
1785 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N;
1786 break;
1787 }
1788 } else {
1789 switch (bs) {
1790 case 1:
1791 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace;
1792 break;
1793 case 2:
1794 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering;
1795 break;
1796 case 3:
1797 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering;
1798 break;
1799 case 4:
1800 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering;
1801 break;
1802 case 5:
1803 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering;
1804 break;
1805 case 6:
1806 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering;
1807 break;
1808 case 7:
1809 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering;
1810 break;
1811 default:
1812 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering;
1813 break;
1814 }
1815 }
1816 PetscFunctionReturn(0);
1817 }
1818
1819 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
1820 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*);
1821
MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat * B)1822 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B)
1823 {
1824 PetscInt n = A->rmap->n;
1825 PetscErrorCode ierr;
1826
1827 PetscFunctionBegin;
1828 #if defined(PETSC_USE_COMPLEX)
1829 if (A->hermitian && !A->symmetric && (ftype == MAT_FACTOR_CHOLESKY||ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported");
1830 #endif
1831
1832 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr);
1833 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr);
1834 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1835 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr);
1836 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
1837
1838 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ;
1839 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ;
1840 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
1841
1842 (*B)->factortype = ftype;
1843 (*B)->useordering = PETSC_TRUE;
1844 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr);
1845 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr);
1846 PetscFunctionReturn(0);
1847 }
1848
1849 /*@C
1850 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored
1851
1852 Not Collective
1853
1854 Input Parameter:
1855 . mat - a MATSEQSBAIJ matrix
1856
1857 Output Parameter:
1858 . array - pointer to the data
1859
1860 Level: intermediate
1861
1862 .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1863 @*/
MatSeqSBAIJGetArray(Mat A,PetscScalar ** array)1864 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array)
1865 {
1866 PetscErrorCode ierr;
1867
1868 PetscFunctionBegin;
1869 ierr = PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
1870 PetscFunctionReturn(0);
1871 }
1872
1873 /*@C
1874 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray()
1875
1876 Not Collective
1877
1878 Input Parameters:
1879 + mat - a MATSEQSBAIJ matrix
1880 - array - pointer to the data
1881
1882 Level: intermediate
1883
1884 .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
1885 @*/
MatSeqSBAIJRestoreArray(Mat A,PetscScalar ** array)1886 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array)
1887 {
1888 PetscErrorCode ierr;
1889
1890 PetscFunctionBegin;
1891 ierr = PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
1892 PetscFunctionReturn(0);
1893 }
1894
1895 /*MC
1896 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices,
1897 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored.
1898
1899 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
1900 can call MatSetOption(Mat, MAT_HERMITIAN).
1901
1902 Options Database Keys:
1903 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions()
1904
1905 Notes:
1906 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not
1907 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use
1908 the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion.
1909
1910 The number of rows in the matrix must be less than or equal to the number of columns
1911
1912 Level: beginner
1913
1914 .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ
1915 M*/
MatCreate_SeqSBAIJ(Mat B)1916 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B)
1917 {
1918 Mat_SeqSBAIJ *b;
1919 PetscErrorCode ierr;
1920 PetscMPIInt size;
1921 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE;
1922
1923 PetscFunctionBegin;
1924 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr);
1925 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
1926
1927 ierr = PetscNewLog(B,&b);CHKERRQ(ierr);
1928 B->data = (void*)b;
1929 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr);
1930
1931 B->ops->destroy = MatDestroy_SeqSBAIJ;
1932 B->ops->view = MatView_SeqSBAIJ;
1933 b->row = NULL;
1934 b->icol = NULL;
1935 b->reallocs = 0;
1936 b->saved_values = NULL;
1937 b->inode.limit = 5;
1938 b->inode.max_limit = 5;
1939
1940 b->roworiented = PETSC_TRUE;
1941 b->nonew = 0;
1942 b->diag = NULL;
1943 b->solve_work = NULL;
1944 b->mult_work = NULL;
1945 B->spptr = NULL;
1946 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2;
1947 b->keepnonzeropattern = PETSC_FALSE;
1948
1949 b->inew = NULL;
1950 b->jnew = NULL;
1951 b->anew = NULL;
1952 b->a2anew = NULL;
1953 b->permute = PETSC_FALSE;
1954
1955 b->ignore_ltriangular = PETSC_TRUE;
1956
1957 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr);
1958
1959 b->getrow_utriangular = PETSC_FALSE;
1960
1961 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr);
1962
1963 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);CHKERRQ(ierr);
1964 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);CHKERRQ(ierr);
1965 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr);
1966 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr);
1967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr);
1968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr);
1969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr);
1970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr);
1971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr);
1972 #if defined(PETSC_HAVE_ELEMENTAL)
1973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr);
1974 #endif
1975 #if defined(PETSC_HAVE_SCALAPACK)
1976 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);CHKERRQ(ierr);
1977 #endif
1978
1979 B->symmetric = PETSC_TRUE;
1980 B->structurally_symmetric = PETSC_TRUE;
1981 B->symmetric_set = PETSC_TRUE;
1982 B->structurally_symmetric_set = PETSC_TRUE;
1983 B->symmetric_eternal = PETSC_TRUE;
1984 #if defined(PETSC_USE_COMPLEX)
1985 B->hermitian = PETSC_FALSE;
1986 B->hermitian_set = PETSC_FALSE;
1987 #else
1988 B->hermitian = PETSC_TRUE;
1989 B->hermitian_set = PETSC_TRUE;
1990 #endif
1991
1992 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr);
1993
1994 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr);
1995 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr);
1996 if (no_unroll) {
1997 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);
1998 }
1999 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr);
2000 if (no_inode) {
2001 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);
2002 }
2003 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr);
2004 ierr = PetscOptionsEnd();CHKERRQ(ierr);
2005 b->inode.use = (PetscBool)(!(no_unroll || no_inode));
2006 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit;
2007 PetscFunctionReturn(0);
2008 }
2009
2010 /*@C
2011 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block
2012 compressed row) format. For good matrix assembly performance the
2013 user should preallocate the matrix storage by setting the parameter nz
2014 (or the array nnz). By setting these parameters accurately, performance
2015 during matrix assembly can be increased by more than a factor of 50.
2016
2017 Collective on Mat
2018
2019 Input Parameters:
2020 + B - the symmetric matrix
2021 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2022 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2023 . nz - number of block nonzeros per block row (same for all rows)
2024 - nnz - array containing the number of block nonzeros in the upper triangular plus
2025 diagonal portion of each block (possibly different for each block row) or NULL
2026
2027 Options Database Keys:
2028 + -mat_no_unroll - uses code that does not unroll the loops in the
2029 block calculations (much slower)
2030 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in
2031
2032 Level: intermediate
2033
2034 Notes:
2035 Specify the preallocated storage with either nz or nnz (not both).
2036 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2037 allocation. See Users-Manual: ch_mat for details.
2038
2039 You can call MatGetInfo() to get information on how effective the preallocation was;
2040 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
2041 You can also run with the option -info and look for messages with the string
2042 malloc in them to see if additional memory allocation was needed.
2043
2044 If the nnz parameter is given then the nz parameter is ignored
2045
2046
2047 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2048 @*/
MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])2049 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
2050 {
2051 PetscErrorCode ierr;
2052
2053 PetscFunctionBegin;
2054 PetscValidHeaderSpecific(B,MAT_CLASSID,1);
2055 PetscValidType(B,1);
2056 PetscValidLogicalCollectiveInt(B,bs,2);
2057 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr);
2058 PetscFunctionReturn(0);
2059 }
2060
2061 /*@C
2062 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values
2063
2064 Input Parameters:
2065 + B - the matrix
2066 . bs - size of block, the blocks are ALWAYS square.
2067 . i - the indices into j for the start of each local row (starts with zero)
2068 . j - the column indices for each local row (starts with zero) these must be sorted for each row
2069 - v - optional values in the matrix
2070
2071 Level: advanced
2072
2073 Notes:
2074 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs
2075 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
2076 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
2077 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2078 block column and the second index is over columns within a block.
2079
2080 Any entries below the diagonal are ignored
2081
2082 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries
2083 and usually the numerical values as well
2084
2085 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ
2086 @*/
MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[],const PetscScalar v[])2087 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
2088 {
2089 PetscErrorCode ierr;
2090
2091 PetscFunctionBegin;
2092 PetscValidHeaderSpecific(B,MAT_CLASSID,1);
2093 PetscValidType(B,1);
2094 PetscValidLogicalCollectiveInt(B,bs,2);
2095 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr);
2096 PetscFunctionReturn(0);
2097 }
2098
2099 /*@C
2100 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block
2101 compressed row) format. For good matrix assembly performance the
2102 user should preallocate the matrix storage by setting the parameter nz
2103 (or the array nnz). By setting these parameters accurately, performance
2104 during matrix assembly can be increased by more than a factor of 50.
2105
2106 Collective
2107
2108 Input Parameters:
2109 + comm - MPI communicator, set to PETSC_COMM_SELF
2110 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
2111 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
2112 . m - number of rows, or number of columns
2113 . nz - number of block nonzeros per block row (same for all rows)
2114 - nnz - array containing the number of block nonzeros in the upper triangular plus
2115 diagonal portion of each block (possibly different for each block row) or NULL
2116
2117 Output Parameter:
2118 . A - the symmetric matrix
2119
2120 Options Database Keys:
2121 + -mat_no_unroll - uses code that does not unroll the loops in the
2122 block calculations (much slower)
2123 - -mat_block_size - size of the blocks to use
2124
2125 Level: intermediate
2126
2127 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
2128 MatXXXXSetPreallocation() paradigm instead of this routine directly.
2129 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
2130
2131 Notes:
2132 The number of rows and columns must be divisible by blocksize.
2133 This matrix type does not support complex Hermitian operation.
2134
2135 Specify the preallocated storage with either nz or nnz (not both).
2136 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
2137 allocation. See Users-Manual: ch_mat for details.
2138
2139 If the nnz parameter is given then the nz parameter is ignored
2140
2141 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ()
2142 @*/
MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat * A)2143 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
2144 {
2145 PetscErrorCode ierr;
2146
2147 PetscFunctionBegin;
2148 ierr = MatCreate(comm,A);CHKERRQ(ierr);
2149 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr);
2150 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr);
2151 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr);
2152 PetscFunctionReturn(0);
2153 }
2154
MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat * B)2155 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
2156 {
2157 Mat C;
2158 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data;
2159 PetscErrorCode ierr;
2160 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2;
2161
2162 PetscFunctionBegin;
2163 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");
2164
2165 *B = NULL;
2166 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr);
2167 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr);
2168 ierr = MatSetBlockSizesFromMats(C,A,A);CHKERRQ(ierr);
2169 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr);
2170 c = (Mat_SeqSBAIJ*)C->data;
2171
2172 C->preallocated = PETSC_TRUE;
2173 C->factortype = A->factortype;
2174 c->row = NULL;
2175 c->icol = NULL;
2176 c->saved_values = NULL;
2177 c->keepnonzeropattern = a->keepnonzeropattern;
2178 C->assembled = PETSC_TRUE;
2179
2180 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr);
2181 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr);
2182 c->bs2 = a->bs2;
2183 c->mbs = a->mbs;
2184 c->nbs = a->nbs;
2185
2186 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2187 c->imax = a->imax;
2188 c->ilen = a->ilen;
2189 c->free_imax_ilen = PETSC_FALSE;
2190 } else {
2191 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr);
2192 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr);
2193 for (i=0; i<mbs; i++) {
2194 c->imax[i] = a->imax[i];
2195 c->ilen[i] = a->ilen[i];
2196 }
2197 c->free_imax_ilen = PETSC_TRUE;
2198 }
2199
2200 /* allocate the matrix space */
2201 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2202 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr);
2203 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr);
2204 c->i = a->i;
2205 c->j = a->j;
2206 c->singlemalloc = PETSC_FALSE;
2207 c->free_a = PETSC_TRUE;
2208 c->free_ij = PETSC_FALSE;
2209 c->parent = A;
2210 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
2211 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
2212 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
2213 } else {
2214 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr);
2215 ierr = PetscArraycpy(c->i,a->i,mbs+1);CHKERRQ(ierr);
2216 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr);
2217 c->singlemalloc = PETSC_TRUE;
2218 c->free_a = PETSC_TRUE;
2219 c->free_ij = PETSC_TRUE;
2220 }
2221 if (mbs > 0) {
2222 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) {
2223 ierr = PetscArraycpy(c->j,a->j,nz);CHKERRQ(ierr);
2224 }
2225 if (cpvalues == MAT_COPY_VALUES) {
2226 ierr = PetscArraycpy(c->a,a->a,bs2*nz);CHKERRQ(ierr);
2227 } else {
2228 ierr = PetscArrayzero(c->a,bs2*nz);CHKERRQ(ierr);
2229 }
2230 if (a->jshort) {
2231 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */
2232 /* if the parent matrix is reassembled, this child matrix will never notice */
2233 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr);
2234 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr);
2235 ierr = PetscArraycpy(c->jshort,a->jshort,nz);CHKERRQ(ierr);
2236
2237 c->free_jshort = PETSC_TRUE;
2238 }
2239 }
2240
2241 c->roworiented = a->roworiented;
2242 c->nonew = a->nonew;
2243
2244 if (a->diag) {
2245 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
2246 c->diag = a->diag;
2247 c->free_diag = PETSC_FALSE;
2248 } else {
2249 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr);
2250 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr);
2251 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
2252 c->free_diag = PETSC_TRUE;
2253 }
2254 }
2255 c->nz = a->nz;
2256 c->maxnz = a->nz; /* Since we allocate exactly the right amount */
2257 c->solve_work = NULL;
2258 c->mult_work = NULL;
2259
2260 *B = C;
2261 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr);
2262 PetscFunctionReturn(0);
2263 }
2264
2265 /* Used for both SeqBAIJ and SeqSBAIJ matrices */
2266 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary
2267
MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)2268 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer)
2269 {
2270 PetscErrorCode ierr;
2271 PetscBool isbinary;
2272
2273 PetscFunctionBegin;
2274 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
2275 if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)mat)->type_name);
2276 ierr = MatLoad_SeqSBAIJ_Binary(mat,viewer);CHKERRQ(ierr);
2277 PetscFunctionReturn(0);
2278 }
2279
2280 /*@
2281 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements
2282 (upper triangular entries in CSR format) provided by the user.
2283
2284 Collective
2285
2286 Input Parameters:
2287 + comm - must be an MPI communicator of size 1
2288 . bs - size of block
2289 . m - number of rows
2290 . n - number of columns
2291 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix
2292 . j - column indices
2293 - a - matrix values
2294
2295 Output Parameter:
2296 . mat - the matrix
2297
2298 Level: advanced
2299
2300 Notes:
2301 The i, j, and a arrays are not copied by this routine, the user must free these arrays
2302 once the matrix is destroyed
2303
2304 You cannot set new nonzero locations into this matrix, that will generate an error.
2305
2306 The i and j indices are 0 based
2307
2308 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1
2309 it is the regular CSR format excluding the lower triangular elements.
2310
2311 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ()
2312
2313 @*/
MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat * mat)2314 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
2315 {
2316 PetscErrorCode ierr;
2317 PetscInt ii;
2318 Mat_SeqSBAIJ *sbaij;
2319
2320 PetscFunctionBegin;
2321 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
2322 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
2323
2324 ierr = MatCreate(comm,mat);CHKERRQ(ierr);
2325 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr);
2326 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr);
2327 ierr = MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr);
2328 sbaij = (Mat_SeqSBAIJ*)(*mat)->data;
2329 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr);
2330 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr);
2331
2332 sbaij->i = i;
2333 sbaij->j = j;
2334 sbaij->a = a;
2335
2336 sbaij->singlemalloc = PETSC_FALSE;
2337 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
2338 sbaij->free_a = PETSC_FALSE;
2339 sbaij->free_ij = PETSC_FALSE;
2340 sbaij->free_imax_ilen = PETSC_TRUE;
2341
2342 for (ii=0; ii<m; ii++) {
2343 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii];
2344 if (PetscUnlikelyDebug(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]);
2345 }
2346 if (PetscDefined(USE_DEBUG)) {
2347 for (ii=0; ii<sbaij->i[m]; ii++) {
2348 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
2349 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]);
2350 }
2351 }
2352
2353 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2354 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2355 PetscFunctionReturn(0);
2356 }
2357
MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat * outmat)2358 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2359 {
2360 PetscErrorCode ierr;
2361 PetscMPIInt size;
2362
2363 PetscFunctionBegin;
2364 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
2365 if (size == 1 && scall == MAT_REUSE_MATRIX) {
2366 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
2367 } else {
2368 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr);
2369 }
2370 PetscFunctionReturn(0);
2371 }
2372