1 /*$Id: baij.c,v 1.202 2000/04/09 04:36:19 bsmith Exp bsmith $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "sys.h" 8 #include "src/mat/impls/baij/seq/baij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 12 /* UGLY, ugly, ugly 13 When MatScalar == Scalar the function MatSetValuesBlocked_SeqBAIJ_MatScalar() does 14 not exist. Otherwise ..._MatScalar() takes matrix elements in single precision and 15 inserts them into the single precision data structure. The function MatSetValuesBlocked_SeqBAIJ() 16 converts the entries into single precision and then calls ..._MatScalar() to put them 17 into the single precision data structures. 18 */ 19 #if defined(PETSC_USE_MAT_SINGLE) 20 extern int MatSetValuesBlocked_SeqBAIJ_MatScalar(Mat,int,int*,int,int*,MatScalar*,InsertMode); 21 #else 22 #define MatSetValuesBlocked_SeqBAIJ_MatScalar MatSetValuesBlocked_SeqBAIJ 23 #endif 24 25 #define CHUNKSIZE 10 26 27 /* 28 Checks for missing diagonals 29 */ 30 #undef __FUNC__ 31 #define __FUNC__ /*<a name=""></a>*/"MatMissingDiagonal_SeqBAIJ" 32 int MatMissingDiagonal_SeqBAIJ(Mat A) 33 { 34 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 35 int *diag,*jj = a->j,i,ierr; 36 37 PetscFunctionBegin; 38 ierr = MatMarkDiagonal_SeqBAIJ(A);CHKERRQ(ierr); 39 diag = a->diag; 40 for (i=0; i<a->mbs; i++) { 41 if (jj[diag[i]] != i) { 42 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 43 } 44 } 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNC__ 49 #define __FUNC__ /*<a name=""></a>*/"MatMarkDiagonal_SeqBAIJ" 50 int MatMarkDiagonal_SeqBAIJ(Mat A) 51 { 52 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 53 int i,j,*diag,m = a->mbs; 54 55 PetscFunctionBegin; 56 if (a->diag) PetscFunctionReturn(0); 57 58 diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(diag); 59 PLogObjectMemory(A,(m+1)*sizeof(int)); 60 for (i=0; i<m; i++) { 61 diag[i] = a->i[i+1]; 62 for (j=a->i[i]; j<a->i[i+1]; j++) { 63 if (a->j[j] == i) { 64 diag[i] = j; 65 break; 66 } 67 } 68 } 69 a->diag = diag; 70 PetscFunctionReturn(0); 71 } 72 73 74 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 75 76 #undef __FUNC__ 77 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ_SeqBAIJ" 78 static int MatGetRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 79 { 80 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 81 int ierr,n = a->mbs,i; 82 83 PetscFunctionBegin; 84 *nn = n; 85 if (!ia) PetscFunctionReturn(0); 86 if (symmetric) { 87 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 88 } else if (oshift == 1) { 89 /* temporarily add 1 to i and j indices */ 90 int nz = a->i[n] + 1; 91 for (i=0; i<nz; i++) a->j[i]++; 92 for (i=0; i<n+1; i++) a->i[i]++; 93 *ia = a->i; *ja = a->j; 94 } else { 95 *ia = a->i; *ja = a->j; 96 } 97 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNC__ 102 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ_SeqBAIJ" 103 static int MatRestoreRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 104 PetscTruth *done) 105 { 106 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 107 int i,n = a->mbs,ierr; 108 109 PetscFunctionBegin; 110 if (!ia) PetscFunctionReturn(0); 111 if (symmetric) { 112 ierr = PetscFree(*ia);CHKERRQ(ierr); 113 ierr = PetscFree(*ja);CHKERRQ(ierr); 114 } else if (oshift == 1) { 115 int nz = a->i[n]; 116 for (i=0; i<nz; i++) a->j[i]--; 117 for (i=0; i<n+1; i++) a->i[i]--; 118 } 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNC__ 123 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize_SeqBAIJ" 124 int MatGetBlockSize_SeqBAIJ(Mat mat,int *bs) 125 { 126 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ*)mat->data; 127 128 PetscFunctionBegin; 129 *bs = baij->bs; 130 PetscFunctionReturn(0); 131 } 132 133 134 #undef __FUNC__ 135 #define __FUNC__ /*<a name=""></a>*/"MatDestroy_SeqBAIJ" 136 int MatDestroy_SeqBAIJ(Mat A) 137 { 138 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 139 int ierr; 140 141 PetscFunctionBegin; 142 if (A->mapping) { 143 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 144 } 145 if (A->bmapping) { 146 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 147 } 148 if (A->rmap) { 149 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 150 } 151 if (A->cmap) { 152 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 153 } 154 #if defined(PETSC_USE_LOG) 155 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 156 #endif 157 ierr = PetscFree(a->a);CHKERRQ(ierr); 158 if (!a->singlemalloc) { 159 ierr = PetscFree(a->i);CHKERRQ(ierr); 160 ierr = PetscFree(a->j);CHKERRQ(ierr); 161 } 162 if (a->row) { 163 ierr = ISDestroy(a->row);CHKERRQ(ierr); 164 } 165 if (a->col) { 166 ierr = ISDestroy(a->col);CHKERRQ(ierr); 167 } 168 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 169 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 170 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 171 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 172 if (a->mult_work) {ierr = PetscFree(a->mult_work);CHKERRQ(ierr);} 173 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 174 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 175 ierr = PetscFree(a);CHKERRQ(ierr); 176 PLogObjectDestroy(A); 177 PetscHeaderDestroy(A); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNC__ 182 #define __FUNC__ /*<a name=""></a>*/"MatSetOption_SeqBAIJ" 183 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 184 { 185 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 186 187 PetscFunctionBegin; 188 if (op == MAT_ROW_ORIENTED) a->roworiented = PETSC_TRUE; 189 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = PETSC_FALSE; 190 else if (op == MAT_COLUMNS_SORTED) a->sorted = PETSC_TRUE; 191 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = PETSC_FALSE; 192 else if (op == MAT_KEEP_ZEROED_ROWS) a->keepzeroedrows = PETSC_TRUE; 193 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 194 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 195 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 196 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 197 else if (op == MAT_ROWS_SORTED || 198 op == MAT_ROWS_UNSORTED || 199 op == MAT_SYMMETRIC || 200 op == MAT_STRUCTURALLY_SYMMETRIC || 201 op == MAT_YES_NEW_DIAGONALS || 202 op == MAT_IGNORE_OFF_PROC_ENTRIES || 203 op == MAT_USE_HASH_TABLE) { 204 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 205 } else if (op == MAT_NO_NEW_DIAGONALS) { 206 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 207 } else { 208 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 209 } 210 PetscFunctionReturn(0); 211 } 212 213 214 #undef __FUNC__ 215 #define __FUNC__ /*<a name=""></a>*/"MatGetSize_SeqBAIJ" 216 int MatGetSize_SeqBAIJ(Mat A,int *m,int *n) 217 { 218 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 219 220 PetscFunctionBegin; 221 if (m) *m = a->m; 222 if (n) *n = a->n; 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNC__ 227 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange_SeqBAIJ" 228 int MatGetOwnershipRange_SeqBAIJ(Mat A,int *m,int *n) 229 { 230 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 231 232 PetscFunctionBegin; 233 if (m) *m = 0; 234 if (n) *n = a->m; 235 PetscFunctionReturn(0); 236 } 237 238 #undef __FUNC__ 239 #define __FUNC__ /*<a name=""></a>*/"MatGetRow_SeqBAIJ" 240 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 241 { 242 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 243 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2; 244 MatScalar *aa,*aa_i; 245 Scalar *v_i; 246 247 PetscFunctionBegin; 248 bs = a->bs; 249 ai = a->i; 250 aj = a->j; 251 aa = a->a; 252 bs2 = a->bs2; 253 254 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 255 256 bn = row/bs; /* Block number */ 257 bp = row % bs; /* Block Position */ 258 M = ai[bn+1] - ai[bn]; 259 *nz = bs*M; 260 261 if (v) { 262 *v = 0; 263 if (*nz) { 264 *v = (Scalar*)PetscMalloc((*nz)*sizeof(Scalar));CHKPTRQ(*v); 265 for (i=0; i<M; i++) { /* for each block in the block row */ 266 v_i = *v + i*bs; 267 aa_i = aa + bs2*(ai[bn] + i); 268 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 269 } 270 } 271 } 272 273 if (idx) { 274 *idx = 0; 275 if (*nz) { 276 *idx = (int*)PetscMalloc((*nz)*sizeof(int));CHKPTRQ(*idx); 277 for (i=0; i<M; i++) { /* for each block in the block row */ 278 idx_i = *idx + i*bs; 279 itmp = bs*aj[ai[bn] + i]; 280 for (j=0; j<bs; j++) {idx_i[j] = itmp++;} 281 } 282 } 283 } 284 PetscFunctionReturn(0); 285 } 286 287 #undef __FUNC__ 288 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRow_SeqBAIJ" 289 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 290 { 291 int ierr; 292 293 PetscFunctionBegin; 294 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 295 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 296 PetscFunctionReturn(0); 297 } 298 299 #undef __FUNC__ 300 #define __FUNC__ /*<a name=""></a>*/"MatTranspose_SeqBAIJ" 301 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 302 { 303 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 304 Mat C; 305 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 306 int *rows,*cols,bs2=a->bs2,refct; 307 Scalar *array; 308 309 PetscFunctionBegin; 310 if (!B && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 311 col = (int*)PetscMalloc((1+nbs)*sizeof(int));CHKPTRQ(col); 312 ierr = PetscMemzero(col,(1+nbs)*sizeof(int));CHKERRQ(ierr); 313 314 #if defined(PETSC_USE_MAT_SINGLE) 315 array = (Scalar *)PetscMalloc(a->bs2*a->nz*sizeof(Scalar));CHKPTRQ(array); 316 for (i=0; i<a->bs2*a->nz; i++) array[i] = (Scalar)a->a[i]; 317 #else 318 array = a->a; 319 #endif 320 321 for (i=0; i<ai[mbs]; i++) col[aj[i]] += 1; 322 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C);CHKERRQ(ierr); 323 ierr = PetscFree(col);CHKERRQ(ierr); 324 rows = (int*)PetscMalloc(2*bs*sizeof(int));CHKPTRQ(rows); 325 cols = rows + bs; 326 for (i=0; i<mbs; i++) { 327 cols[0] = i*bs; 328 for (k=1; k<bs; k++) cols[k] = cols[k-1] + 1; 329 len = ai[i+1] - ai[i]; 330 for (j=0; j<len; j++) { 331 rows[0] = (*aj++)*bs; 332 for (k=1; k<bs; k++) rows[k] = rows[k-1] + 1; 333 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES);CHKERRQ(ierr); 334 array += bs2; 335 } 336 } 337 ierr = PetscFree(rows);CHKERRQ(ierr); 338 #if defined(PETSC_USE_MAT_SINGLE) 339 ierr = PetscFree(array); 340 #endif 341 342 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 343 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 344 345 if (B) { 346 *B = C; 347 } else { 348 PetscOps *Abops; 349 MatOps Aops; 350 351 /* This isn't really an in-place transpose */ 352 ierr = PetscFree(a->a);CHKERRQ(ierr); 353 if (!a->singlemalloc) { 354 ierr = PetscFree(a->i);CHKERRQ(ierr); 355 ierr = PetscFree(a->j);CHKERRQ(ierr); 356 } 357 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 358 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 359 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 360 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 361 ierr = PetscFree(a);CHKERRQ(ierr); 362 363 364 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 365 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 366 367 /* 368 This is horrible, horrible code. We need to keep the 369 A pointers for the bops and ops but copy everything 370 else from C. 371 */ 372 Abops = A->bops; 373 Aops = A->ops; 374 refct = A->refct; 375 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 376 A->bops = Abops; 377 A->ops = Aops; 378 A->qlist = 0; 379 A->refct = refct; 380 381 PetscHeaderDestroy(C); 382 } 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNC__ 387 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqBAIJ_Binary" 388 static int MatView_SeqBAIJ_Binary(Mat A,Viewer viewer) 389 { 390 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 391 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 392 Scalar *aa; 393 FILE *file; 394 395 PetscFunctionBegin; 396 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 397 col_lens = (int*)PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 398 col_lens[0] = MAT_COOKIE; 399 400 col_lens[1] = a->m; 401 col_lens[2] = a->n; 402 col_lens[3] = a->nz*bs2; 403 404 /* store lengths of each row and write (including header) to file */ 405 count = 0; 406 for (i=0; i<a->mbs; i++) { 407 for (j=0; j<bs; j++) { 408 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 409 } 410 } 411 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1);CHKERRQ(ierr); 412 ierr = PetscFree(col_lens);CHKERRQ(ierr); 413 414 /* store column indices (zero start index) */ 415 jj = (int*)PetscMalloc((a->nz+1)*bs2*sizeof(int));CHKPTRQ(jj); 416 count = 0; 417 for (i=0; i<a->mbs; i++) { 418 for (j=0; j<bs; j++) { 419 for (k=a->i[i]; k<a->i[i+1]; k++) { 420 for (l=0; l<bs; l++) { 421 jj[count++] = bs*a->j[k] + l; 422 } 423 } 424 } 425 } 426 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0);CHKERRQ(ierr); 427 ierr = PetscFree(jj);CHKERRQ(ierr); 428 429 /* store nonzero values */ 430 aa = (Scalar*)PetscMalloc((a->nz+1)*bs2*sizeof(Scalar));CHKPTRQ(aa); 431 count = 0; 432 for (i=0; i<a->mbs; i++) { 433 for (j=0; j<bs; j++) { 434 for (k=a->i[i]; k<a->i[i+1]; k++) { 435 for (l=0; l<bs; l++) { 436 aa[count++] = a->a[bs2*k + l*bs + j]; 437 } 438 } 439 } 440 } 441 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 442 ierr = PetscFree(aa);CHKERRQ(ierr); 443 444 ierr = ViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 445 if (file) { 446 fprintf(file,"-matload_block_size %d\n",a->bs); 447 } 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNC__ 452 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqBAIJ_ASCII" 453 static int MatView_SeqBAIJ_ASCII(Mat A,Viewer viewer) 454 { 455 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 456 int ierr,i,j,format,bs = a->bs,k,l,bs2=a->bs2; 457 char *outputname; 458 459 PetscFunctionBegin; 460 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 461 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 462 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 463 ierr = ViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 464 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 465 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 466 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 467 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 468 for (i=0; i<a->mbs; i++) { 469 for (j=0; j<bs; j++) { 470 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 471 for (k=a->i[i]; k<a->i[i+1]; k++) { 472 for (l=0; l<bs; l++) { 473 #if defined(PETSC_USE_COMPLEX) 474 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 475 ierr = ViewerASCIIPrintf(viewer," %d %g + %gi",bs*a->j[k]+l, 476 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 477 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 478 ierr = ViewerASCIIPrintf(viewer," %d %g - %gi",bs*a->j[k]+l, 479 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 480 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 481 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 482 } 483 #else 484 if (a->a[bs2*k + l*bs + j] != 0.0) { 485 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 486 } 487 #endif 488 } 489 } 490 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 491 } 492 } 493 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 494 } else { 495 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 496 for (i=0; i<a->mbs; i++) { 497 for (j=0; j<bs; j++) { 498 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 499 for (k=a->i[i]; k<a->i[i+1]; k++) { 500 for (l=0; l<bs; l++) { 501 #if defined(PETSC_USE_COMPLEX) 502 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 503 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 504 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 505 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 506 ierr = ViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 507 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 508 } else { 509 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 510 } 511 #else 512 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 513 #endif 514 } 515 } 516 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 517 } 518 } 519 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 520 } 521 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 522 PetscFunctionReturn(0); 523 } 524 525 #undef __FUNC__ 526 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqBAIJ_Draw_Zoom" 527 static int MatView_SeqBAIJ_Draw_Zoom(Draw draw,void *Aa) 528 { 529 Mat A = (Mat) Aa; 530 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 531 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2; 532 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 533 MatScalar *aa; 534 Viewer viewer; 535 536 PetscFunctionBegin; 537 538 /* still need to add support for contour plot of nonzeros; see MatView_SeqAIJ_Draw_Zoom()*/ 539 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 540 541 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 542 543 /* loop over matrix elements drawing boxes */ 544 color = DRAW_BLUE; 545 for (i=0,row=0; i<mbs; i++,row+=bs) { 546 for (j=a->i[i]; j<a->i[i+1]; j++) { 547 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 548 x_l = a->j[j]*bs; x_r = x_l + 1.0; 549 aa = a->a + j*bs2; 550 for (k=0; k<bs; k++) { 551 for (l=0; l<bs; l++) { 552 if (PetscRealPart(*aa++) >= 0.) continue; 553 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 554 } 555 } 556 } 557 } 558 color = DRAW_CYAN; 559 for (i=0,row=0; i<mbs; i++,row+=bs) { 560 for (j=a->i[i]; j<a->i[i+1]; j++) { 561 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 562 x_l = a->j[j]*bs; x_r = x_l + 1.0; 563 aa = a->a + j*bs2; 564 for (k=0; k<bs; k++) { 565 for (l=0; l<bs; l++) { 566 if (PetscRealPart(*aa++) != 0.) continue; 567 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 568 } 569 } 570 } 571 } 572 573 color = DRAW_RED; 574 for (i=0,row=0; i<mbs; i++,row+=bs) { 575 for (j=a->i[i]; j<a->i[i+1]; j++) { 576 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 577 x_l = a->j[j]*bs; x_r = x_l + 1.0; 578 aa = a->a + j*bs2; 579 for (k=0; k<bs; k++) { 580 for (l=0; l<bs; l++) { 581 if (PetscRealPart(*aa++) <= 0.) continue; 582 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 583 } 584 } 585 } 586 } 587 PetscFunctionReturn(0); 588 } 589 590 #undef __FUNC__ 591 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqBAIJ_Draw" 592 static int MatView_SeqBAIJ_Draw(Mat A,Viewer viewer) 593 { 594 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 595 int ierr; 596 PetscReal xl,yl,xr,yr,w,h; 597 Draw draw; 598 PetscTruth isnull; 599 600 PetscFunctionBegin; 601 602 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 603 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 604 605 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 606 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 607 xr += w; yr += h; xl = -w; yl = -h; 608 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 609 ierr = DrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 610 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 611 PetscFunctionReturn(0); 612 } 613 614 #undef __FUNC__ 615 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqBAIJ" 616 int MatView_SeqBAIJ(Mat A,Viewer viewer) 617 { 618 int ierr; 619 PetscTruth issocket,isascii,isbinary,isdraw; 620 621 PetscFunctionBegin; 622 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 623 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 624 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 625 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 626 if (issocket) { 627 SETERRQ(PETSC_ERR_SUP,0,"Socket viewer not supported"); 628 } else if (isascii){ 629 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 630 } else if (isbinary) { 631 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 632 } else if (isdraw) { 633 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 634 } else { 635 SETERRQ1(1,1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 636 } 637 PetscFunctionReturn(0); 638 } 639 640 641 #undef __FUNC__ 642 #define __FUNC__ /*<a name=""></a>*/"MatGetValues_SeqBAIJ" 643 int MatGetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 644 { 645 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 646 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 647 int *ai = a->i,*ailen = a->ilen; 648 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 649 MatScalar *ap,*aa = a->a,zero = 0.0; 650 651 PetscFunctionBegin; 652 for (k=0; k<m; k++) { /* loop over rows */ 653 row = im[k]; brow = row/bs; 654 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 655 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 656 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 657 nrow = ailen[brow]; 658 for (l=0; l<n; l++) { /* loop over columns */ 659 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 660 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 661 col = in[l] ; 662 bcol = col/bs; 663 cidx = col%bs; 664 ridx = row%bs; 665 high = nrow; 666 low = 0; /* assume unsorted */ 667 while (high-low > 5) { 668 t = (low+high)/2; 669 if (rp[t] > bcol) high = t; 670 else low = t; 671 } 672 for (i=low; i<high; i++) { 673 if (rp[i] > bcol) break; 674 if (rp[i] == bcol) { 675 *v++ = ap[bs2*i+bs*cidx+ridx]; 676 goto finished; 677 } 678 } 679 *v++ = zero; 680 finished:; 681 } 682 } 683 PetscFunctionReturn(0); 684 } 685 686 #if defined(PETSC_USE_MAT_SINGLE) 687 #undef __FUNC__ 688 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlocked_SeqBAIJ" 689 int MatSetValuesBlocked_SeqBAIJ(Mat mat,int m,int *im,int n,int *in,Scalar *v,InsertMode addv) 690 { 691 int ierr,i,N = m*n*((Mat_SeqBAIJ*)mat->data)->bs2; 692 MatScalar *vsingle = (MatScalar*)v; 693 694 PetscFunctionBegin; 695 for (i=0; i<N; i++) { 696 vsingle[i] = v[i]; 697 } 698 ierr = MatSetValuesBlocked_SeqBAIJ_MatScalar(mat,m,im,n,in,vsingle,addv);CHKERRQ(ierr); 699 PetscFunctionReturn(0); 700 } 701 #endif 702 703 704 #undef __FUNC__ 705 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlocked_SeqBAIJ" 706 int MatSetValuesBlocked_SeqBAIJ_MatScalar(Mat A,int m,int *im,int n,int *in,MatScalar *v,InsertMode is) 707 { 708 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 709 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 710 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 711 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval,ierr; 712 MatScalar *value = v,*ap,*aa = a->a,*bap; 713 714 PetscFunctionBegin; 715 if (roworiented) { 716 stepval = (n-1)*bs; 717 } else { 718 stepval = (m-1)*bs; 719 } 720 for (k=0; k<m; k++) { /* loop over added rows */ 721 row = im[k]; 722 if (row < 0) continue; 723 #if defined(PETSC_USE_BOPT_g) 724 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 725 #endif 726 rp = aj + ai[row]; 727 ap = aa + bs2*ai[row]; 728 rmax = imax[row]; 729 nrow = ailen[row]; 730 low = 0; 731 for (l=0; l<n; l++) { /* loop over added columns */ 732 if (in[l] < 0) continue; 733 #if defined(PETSC_USE_BOPT_g) 734 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 735 #endif 736 col = in[l]; 737 if (roworiented) { 738 value = v + k*(stepval+bs)*bs + l*bs; 739 } else { 740 value = v + l*(stepval+bs)*bs + k*bs; 741 } 742 if (!sorted) low = 0; high = nrow; 743 while (high-low > 7) { 744 t = (low+high)/2; 745 if (rp[t] > col) high = t; 746 else low = t; 747 } 748 for (i=low; i<high; i++) { 749 if (rp[i] > col) break; 750 if (rp[i] == col) { 751 bap = ap + bs2*i; 752 if (roworiented) { 753 if (is == ADD_VALUES) { 754 for (ii=0; ii<bs; ii++,value+=stepval) { 755 for (jj=ii; jj<bs2; jj+=bs) { 756 bap[jj] += *value++; 757 } 758 } 759 } else { 760 for (ii=0; ii<bs; ii++,value+=stepval) { 761 for (jj=ii; jj<bs2; jj+=bs) { 762 bap[jj] = *value++; 763 } 764 } 765 } 766 } else { 767 if (is == ADD_VALUES) { 768 for (ii=0; ii<bs; ii++,value+=stepval) { 769 for (jj=0; jj<bs; jj++) { 770 *bap++ += *value++; 771 } 772 } 773 } else { 774 for (ii=0; ii<bs; ii++,value+=stepval) { 775 for (jj=0; jj<bs; jj++) { 776 *bap++ = *value++; 777 } 778 } 779 } 780 } 781 goto noinsert2; 782 } 783 } 784 if (nonew == 1) goto noinsert2; 785 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 786 if (nrow >= rmax) { 787 /* there is no extra room in row, therefore enlarge */ 788 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 789 MatScalar *new_a; 790 791 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 792 793 /* malloc new storage space */ 794 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 795 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 796 new_j = (int*)(new_a + bs2*new_nz); 797 new_i = new_j + new_nz; 798 799 /* copy over old data into new slots */ 800 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 801 for (ii=row+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 802 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 803 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 804 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 805 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 806 ierr = PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 807 ierr = PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 808 /* free up old matrix storage */ 809 ierr = PetscFree(a->a);CHKERRQ(ierr); 810 if (!a->singlemalloc) { 811 ierr = PetscFree(a->i);CHKERRQ(ierr); 812 ierr = PetscFree(a->j);CHKERRQ(ierr); 813 } 814 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 815 a->singlemalloc = PETSC_TRUE; 816 817 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 818 rmax = imax[row] = imax[row] + CHUNKSIZE; 819 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 820 a->maxnz += bs2*CHUNKSIZE; 821 a->reallocs++; 822 a->nz++; 823 } 824 N = nrow++ - 1; 825 /* shift up all the later entries in this row */ 826 for (ii=N; ii>=i; ii--) { 827 rp[ii+1] = rp[ii]; 828 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 829 } 830 if (N >= i) { 831 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 832 } 833 rp[i] = col; 834 bap = ap + bs2*i; 835 if (roworiented) { 836 for (ii=0; ii<bs; ii++,value+=stepval) { 837 for (jj=ii; jj<bs2; jj+=bs) { 838 bap[jj] = *value++; 839 } 840 } 841 } else { 842 for (ii=0; ii<bs; ii++,value+=stepval) { 843 for (jj=0; jj<bs; jj++) { 844 *bap++ = *value++; 845 } 846 } 847 } 848 noinsert2:; 849 low = i; 850 } 851 ailen[row] = nrow; 852 } 853 PetscFunctionReturn(0); 854 } 855 856 857 #undef __FUNC__ 858 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyEnd_SeqBAIJ" 859 int MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode) 860 { 861 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 862 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 863 int m = a->m,*ip,N,*ailen = a->ilen; 864 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 865 MatScalar *aa = a->a,*ap; 866 867 PetscFunctionBegin; 868 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 869 870 if (m) rmax = ailen[0]; 871 for (i=1; i<mbs; i++) { 872 /* move each row back by the amount of empty slots (fshift) before it*/ 873 fshift += imax[i-1] - ailen[i-1]; 874 rmax = PetscMax(rmax,ailen[i]); 875 if (fshift) { 876 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 877 N = ailen[i]; 878 for (j=0; j<N; j++) { 879 ip[j-fshift] = ip[j]; 880 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 881 } 882 } 883 ai[i] = ai[i-1] + ailen[i-1]; 884 } 885 if (mbs) { 886 fshift += imax[mbs-1] - ailen[mbs-1]; 887 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 888 } 889 /* reset ilen and imax for each row */ 890 for (i=0; i<mbs; i++) { 891 ailen[i] = imax[i] = ai[i+1] - ai[i]; 892 } 893 a->nz = ai[mbs]; 894 895 /* diagonals may have moved, so kill the diagonal pointers */ 896 if (fshift && a->diag) { 897 ierr = PetscFree(a->diag);CHKERRQ(ierr); 898 PLogObjectMemory(A,-(m+1)*sizeof(int)); 899 a->diag = 0; 900 } 901 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 902 m,a->n,a->bs,fshift*bs2,a->nz*bs2); 903 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Number of mallocs during MatSetValues is %d\n", 904 a->reallocs); 905 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 906 a->reallocs = 0; 907 A->info.nz_unneeded = (PetscReal)fshift*bs2; 908 909 PetscFunctionReturn(0); 910 } 911 912 913 914 /* 915 This function returns an array of flags which indicate the locations of contiguous 916 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 917 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 918 Assume: sizes should be long enough to hold all the values. 919 */ 920 #undef __FUNC__ 921 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows_SeqBAIJ_Check_Blocks" 922 static int MatZeroRows_SeqBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 923 { 924 int i,j,k,row; 925 PetscTruth flg; 926 927 PetscFunctionBegin; 928 for (i=0,j=0; i<n; j++) { 929 row = idx[i]; 930 if (row%bs!=0) { /* Not the begining of a block */ 931 sizes[j] = 1; 932 i++; 933 } else if (i+bs > n) { /* complete block doesn't exist (at idx end) */ 934 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 935 i++; 936 } else { /* Begining of the block, so check if the complete block exists */ 937 flg = PETSC_TRUE; 938 for (k=1; k<bs; k++) { 939 if (row+k != idx[i+k]) { /* break in the block */ 940 flg = PETSC_FALSE; 941 break; 942 } 943 } 944 if (flg == PETSC_TRUE) { /* No break in the bs */ 945 sizes[j] = bs; 946 i+= bs; 947 } else { 948 sizes[j] = 1; 949 i++; 950 } 951 } 952 } 953 *bs_max = j; 954 PetscFunctionReturn(0); 955 } 956 957 #undef __FUNC__ 958 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows_SeqBAIJ" 959 int MatZeroRows_SeqBAIJ(Mat A,IS is,Scalar *diag) 960 { 961 Mat_SeqBAIJ *baij=(Mat_SeqBAIJ*)A->data; 962 int ierr,i,j,k,count,is_n,*is_idx,*rows; 963 int bs=baij->bs,bs2=baij->bs2,*sizes,row,bs_max; 964 Scalar zero = 0.0; 965 MatScalar *aa; 966 967 PetscFunctionBegin; 968 /* Make a copy of the IS and sort it */ 969 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 970 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 971 972 /* allocate memory for rows,sizes */ 973 rows = (int*)PetscMalloc((3*is_n+1)*sizeof(int));CHKPTRQ(rows); 974 sizes = rows + is_n; 975 976 /* initialize copy IS valurs to rows, and sort them */ 977 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 978 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 979 if (baij->keepzeroedrows) { 980 for (i=0; i<is_n; i++) { sizes[i] = 1; } 981 bs_max = is_n; 982 } else { 983 ierr = MatZeroRows_SeqBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 984 } 985 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 986 987 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 988 row = rows[j]; 989 if (row < 0 || row > baij->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"row %d out of range",row); 990 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 991 aa = baij->a + baij->i[row/bs]*bs2 + (row%bs); 992 if (sizes[i] == bs && !baij->keepzeroedrows) { 993 if (diag) { 994 if (baij->ilen[row/bs] > 0) { 995 baij->ilen[row/bs] = 1; 996 baij->j[baij->i[row/bs]] = row/bs; 997 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 998 } 999 /* Now insert all the diagoanl values for this bs */ 1000 for (k=0; k<bs; k++) { 1001 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 1002 } 1003 } else { /* (!diag) */ 1004 baij->ilen[row/bs] = 0; 1005 } /* end (!diag) */ 1006 } else { /* (sizes[i] != bs) */ 1007 #if defined (PETSC_USE_DEBUG) 1008 if (sizes[i] != 1) SETERRQ(1,0,"Internal Error. Value should be 1"); 1009 #endif 1010 for (k=0; k<count; k++) { 1011 aa[0] = zero; 1012 aa+=bs; 1013 } 1014 if (diag) { 1015 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 1016 } 1017 } 1018 } 1019 1020 ierr = PetscFree(rows);CHKERRQ(ierr); 1021 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1022 PetscFunctionReturn(0); 1023 } 1024 1025 #undef __FUNC__ 1026 #define __FUNC__ /*<a name=""></a>*/"MatSetValues_SeqBAIJ" 1027 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 1028 { 1029 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1030 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 1031 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 1032 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 1033 int ridx,cidx,bs2=a->bs2,ierr; 1034 MatScalar *ap,value,*aa=a->a,*bap; 1035 1036 PetscFunctionBegin; 1037 for (k=0; k<m; k++) { /* loop over added rows */ 1038 row = im[k]; brow = row/bs; 1039 if (row < 0) continue; 1040 #if defined(PETSC_USE_BOPT_g) 1041 if (row >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large: row %d max %d",row,a->m); 1042 #endif 1043 rp = aj + ai[brow]; 1044 ap = aa + bs2*ai[brow]; 1045 rmax = imax[brow]; 1046 nrow = ailen[brow]; 1047 low = 0; 1048 for (l=0; l<n; l++) { /* loop over added columns */ 1049 if (in[l] < 0) continue; 1050 #if defined(PETSC_USE_BOPT_g) 1051 if (in[l] >= a->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large: col %d max %d",in[l],a->n); 1052 #endif 1053 col = in[l]; bcol = col/bs; 1054 ridx = row % bs; cidx = col % bs; 1055 if (roworiented) { 1056 value = v[l + k*n]; 1057 } else { 1058 value = v[k + l*m]; 1059 } 1060 if (!sorted) low = 0; high = nrow; 1061 while (high-low > 7) { 1062 t = (low+high)/2; 1063 if (rp[t] > bcol) high = t; 1064 else low = t; 1065 } 1066 for (i=low; i<high; i++) { 1067 if (rp[i] > bcol) break; 1068 if (rp[i] == bcol) { 1069 bap = ap + bs2*i + bs*cidx + ridx; 1070 if (is == ADD_VALUES) *bap += value; 1071 else *bap = value; 1072 goto noinsert1; 1073 } 1074 } 1075 if (nonew == 1) goto noinsert1; 1076 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1077 if (nrow >= rmax) { 1078 /* there is no extra room in row, therefore enlarge */ 1079 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 1080 MatScalar *new_a; 1081 1082 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1083 1084 /* Malloc new storage space */ 1085 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 1086 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 1087 new_j = (int*)(new_a + bs2*new_nz); 1088 new_i = new_j + new_nz; 1089 1090 /* copy over old data into new slots */ 1091 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 1092 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 1093 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 1094 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 1095 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 1096 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1097 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 1098 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 1099 /* free up old matrix storage */ 1100 ierr = PetscFree(a->a);CHKERRQ(ierr); 1101 if (!a->singlemalloc) { 1102 ierr = PetscFree(a->i);CHKERRQ(ierr); 1103 ierr = PetscFree(a->j);CHKERRQ(ierr); 1104 } 1105 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 1106 a->singlemalloc = PETSC_TRUE; 1107 1108 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1109 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1110 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 1111 a->maxnz += bs2*CHUNKSIZE; 1112 a->reallocs++; 1113 a->nz++; 1114 } 1115 N = nrow++ - 1; 1116 /* shift up all the later entries in this row */ 1117 for (ii=N; ii>=i; ii--) { 1118 rp[ii+1] = rp[ii]; 1119 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1120 } 1121 if (N>=i) { 1122 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1123 } 1124 rp[i] = bcol; 1125 ap[bs2*i + bs*cidx + ridx] = value; 1126 noinsert1:; 1127 low = i; 1128 } 1129 ailen[brow] = nrow; 1130 } 1131 PetscFunctionReturn(0); 1132 } 1133 1134 extern int MatLUFactorSymbolic_SeqBAIJ(Mat,IS,IS,PetscReal,Mat*); 1135 extern int MatLUFactor_SeqBAIJ(Mat,IS,IS,PetscReal); 1136 extern int MatIncreaseOverlap_SeqBAIJ(Mat,int,IS*,int); 1137 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 1138 extern int MatGetSubMatrices_SeqBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 1139 extern int MatMultTranspose_SeqBAIJ(Mat,Vec,Vec); 1140 extern int MatMultTransposeAdd_SeqBAIJ(Mat,Vec,Vec,Vec); 1141 extern int MatScale_SeqBAIJ(Scalar*,Mat); 1142 extern int MatNorm_SeqBAIJ(Mat,NormType,PetscReal *); 1143 extern int MatEqual_SeqBAIJ(Mat,Mat,PetscTruth*); 1144 extern int MatGetDiagonal_SeqBAIJ(Mat,Vec); 1145 extern int MatDiagonalScale_SeqBAIJ(Mat,Vec,Vec); 1146 extern int MatGetInfo_SeqBAIJ(Mat,MatInfoType,MatInfo *); 1147 extern int MatZeroEntries_SeqBAIJ(Mat); 1148 1149 extern int MatSolve_SeqBAIJ_N(Mat,Vec,Vec); 1150 extern int MatSolve_SeqBAIJ_1(Mat,Vec,Vec); 1151 extern int MatSolve_SeqBAIJ_2(Mat,Vec,Vec); 1152 extern int MatSolve_SeqBAIJ_3(Mat,Vec,Vec); 1153 extern int MatSolve_SeqBAIJ_4(Mat,Vec,Vec); 1154 extern int MatSolve_SeqBAIJ_5(Mat,Vec,Vec); 1155 extern int MatSolve_SeqBAIJ_6(Mat,Vec,Vec); 1156 extern int MatSolve_SeqBAIJ_7(Mat,Vec,Vec); 1157 extern int MatSolveTranspose_SeqBAIJ_7(Mat,Vec,Vec); 1158 extern int MatSolveTranspose_SeqBAIJ_6(Mat,Vec,Vec); 1159 extern int MatSolveTranspose_SeqBAIJ_5(Mat,Vec,Vec); 1160 extern int MatSolveTranspose_SeqBAIJ_4(Mat,Vec,Vec); 1161 extern int MatSolveTranspose_SeqBAIJ_3(Mat,Vec,Vec); 1162 extern int MatSolveTranspose_SeqBAIJ_2(Mat,Vec,Vec); 1163 extern int MatSolveTranspose_SeqBAIJ_1(Mat,Vec,Vec); 1164 1165 extern int MatLUFactorNumeric_SeqBAIJ_N(Mat,Mat*); 1166 extern int MatLUFactorNumeric_SeqBAIJ_1(Mat,Mat*); 1167 extern int MatLUFactorNumeric_SeqBAIJ_2(Mat,Mat*); 1168 extern int MatLUFactorNumeric_SeqBAIJ_3(Mat,Mat*); 1169 extern int MatLUFactorNumeric_SeqBAIJ_4(Mat,Mat*); 1170 extern int MatLUFactorNumeric_SeqBAIJ_5(Mat,Mat*); 1171 extern int MatLUFactorNumeric_SeqBAIJ_6(Mat,Mat*); 1172 1173 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1174 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1175 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1176 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1177 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1178 extern int MatMult_SeqBAIJ_6(Mat,Vec,Vec); 1179 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1180 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1181 1182 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1183 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1184 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1185 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1186 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1187 extern int MatMultAdd_SeqBAIJ_6(Mat,Vec,Vec,Vec); 1188 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1189 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1190 1191 #undef __FUNC__ 1192 #define __FUNC__ /*<a name=""></a>*/"MatILUFactor_SeqBAIJ" 1193 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1194 { 1195 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)inA->data; 1196 Mat outA; 1197 int ierr; 1198 PetscTruth row_identity,col_identity; 1199 1200 PetscFunctionBegin; 1201 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,0,"Only levels = 0 supported for in-place ILU"); 1202 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1203 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1204 if (!row_identity || !col_identity) { 1205 SETERRQ(1,1,"Row and column permutations must be identity for in-place ILU"); 1206 } 1207 1208 outA = inA; 1209 inA->factor = FACTOR_LU; 1210 1211 if (!a->diag) { 1212 ierr = MatMarkDiagonal_SeqBAIJ(inA);CHKERRQ(ierr); 1213 } 1214 /* 1215 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 1216 for ILU(0) factorization with natural ordering 1217 */ 1218 switch (a->bs) { 1219 case 1: 1220 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1221 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1222 case 2: 1223 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 1224 inA->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 1225 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1226 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1227 break; 1228 case 3: 1229 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 1230 inA->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; 1231 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; 1232 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1233 break; 1234 case 4: 1235 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1236 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1237 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1238 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1239 break; 1240 case 5: 1241 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1242 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1243 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1244 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1245 break; 1246 case 6: 1247 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1248 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1249 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1250 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1251 break; 1252 case 7: 1253 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1254 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1255 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1256 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1257 break; 1258 default: 1259 a->row = row; 1260 a->col = col; 1261 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1262 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1263 1264 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1265 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1266 PLogObjectParent(inA,a->icol); 1267 1268 if (!a->solve_work) { 1269 a->solve_work = (Scalar*)PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1270 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1271 } 1272 } 1273 1274 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1275 1276 PetscFunctionReturn(0); 1277 } 1278 #undef __FUNC__ 1279 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp_SeqBAIJ" 1280 int MatPrintHelp_SeqBAIJ(Mat A) 1281 { 1282 static PetscTruth called = PETSC_FALSE; 1283 MPI_Comm comm = A->comm; 1284 int ierr; 1285 1286 PetscFunctionBegin; 1287 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1288 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1289 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1290 PetscFunctionReturn(0); 1291 } 1292 1293 EXTERN_C_BEGIN 1294 #undef __FUNC__ 1295 #define __FUNC__ /*<a name=""></a>*/"MatSeqBAIJSetColumnIndices_SeqBAIJ" 1296 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1297 { 1298 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1299 int i,nz,n; 1300 1301 PetscFunctionBegin; 1302 nz = baij->maxnz; 1303 n = baij->n; 1304 for (i=0; i<nz; i++) { 1305 baij->j[i] = indices[i]; 1306 } 1307 baij->nz = nz; 1308 for (i=0; i<n; i++) { 1309 baij->ilen[i] = baij->imax[i]; 1310 } 1311 1312 PetscFunctionReturn(0); 1313 } 1314 EXTERN_C_END 1315 1316 #undef __FUNC__ 1317 #define __FUNC__ /*<a name=""></a>*/"MatSeqBAIJSetColumnIndices" 1318 /*@ 1319 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1320 in the matrix. 1321 1322 Input Parameters: 1323 + mat - the SeqBAIJ matrix 1324 - indices - the column indices 1325 1326 Level: advanced 1327 1328 Notes: 1329 This can be called if you have precomputed the nonzero structure of the 1330 matrix and want to provide it to the matrix object to improve the performance 1331 of the MatSetValues() operation. 1332 1333 You MUST have set the correct numbers of nonzeros per row in the call to 1334 MatCreateSeqBAIJ(). 1335 1336 MUST be called before any calls to MatSetValues(); 1337 1338 @*/ 1339 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1340 { 1341 int ierr,(*f)(Mat,int *); 1342 1343 PetscFunctionBegin; 1344 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1345 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1346 if (f) { 1347 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1348 } else { 1349 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1350 } 1351 PetscFunctionReturn(0); 1352 } 1353 1354 /* -------------------------------------------------------------------*/ 1355 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1356 MatGetRow_SeqBAIJ, 1357 MatRestoreRow_SeqBAIJ, 1358 MatMult_SeqBAIJ_N, 1359 MatMultAdd_SeqBAIJ_N, 1360 MatMultTranspose_SeqBAIJ, 1361 MatMultTransposeAdd_SeqBAIJ, 1362 MatSolve_SeqBAIJ_N, 1363 0, 1364 0, 1365 0, 1366 MatLUFactor_SeqBAIJ, 1367 0, 1368 0, 1369 MatTranspose_SeqBAIJ, 1370 MatGetInfo_SeqBAIJ, 1371 MatEqual_SeqBAIJ, 1372 MatGetDiagonal_SeqBAIJ, 1373 MatDiagonalScale_SeqBAIJ, 1374 MatNorm_SeqBAIJ, 1375 0, 1376 MatAssemblyEnd_SeqBAIJ, 1377 0, 1378 MatSetOption_SeqBAIJ, 1379 MatZeroEntries_SeqBAIJ, 1380 MatZeroRows_SeqBAIJ, 1381 MatLUFactorSymbolic_SeqBAIJ, 1382 MatLUFactorNumeric_SeqBAIJ_N, 1383 0, 1384 0, 1385 MatGetSize_SeqBAIJ, 1386 MatGetSize_SeqBAIJ, 1387 MatGetOwnershipRange_SeqBAIJ, 1388 MatILUFactorSymbolic_SeqBAIJ, 1389 0, 1390 0, 1391 0, 1392 MatDuplicate_SeqBAIJ, 1393 0, 1394 0, 1395 MatILUFactor_SeqBAIJ, 1396 0, 1397 0, 1398 MatGetSubMatrices_SeqBAIJ, 1399 MatIncreaseOverlap_SeqBAIJ, 1400 MatGetValues_SeqBAIJ, 1401 0, 1402 MatPrintHelp_SeqBAIJ, 1403 MatScale_SeqBAIJ, 1404 0, 1405 0, 1406 0, 1407 MatGetBlockSize_SeqBAIJ, 1408 MatGetRowIJ_SeqBAIJ, 1409 MatRestoreRowIJ_SeqBAIJ, 1410 0, 1411 0, 1412 0, 1413 0, 1414 0, 1415 0, 1416 MatSetValuesBlocked_SeqBAIJ, 1417 MatGetSubMatrix_SeqBAIJ, 1418 0, 1419 0, 1420 MatGetMaps_Petsc}; 1421 1422 EXTERN_C_BEGIN 1423 #undef __FUNC__ 1424 #define __FUNC__ /*<a name=""></a>*/"MatStoreValues_SeqBAIJ" 1425 int MatStoreValues_SeqBAIJ(Mat mat) 1426 { 1427 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1428 int nz = aij->i[aij->m]*aij->bs*aij->bs2; 1429 int ierr; 1430 1431 PetscFunctionBegin; 1432 if (aij->nonew != 1) { 1433 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1434 } 1435 1436 /* allocate space for values if not already there */ 1437 if (!aij->saved_values) { 1438 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1439 } 1440 1441 /* copy values over */ 1442 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1443 PetscFunctionReturn(0); 1444 } 1445 EXTERN_C_END 1446 1447 EXTERN_C_BEGIN 1448 #undef __FUNC__ 1449 #define __FUNC__ /*<a name=""></a>*/"MatRetrieveValues_SeqBAIJ" 1450 int MatRetrieveValues_SeqBAIJ(Mat mat) 1451 { 1452 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1453 int nz = aij->i[aij->m]*aij->bs*aij->bs2,ierr; 1454 1455 PetscFunctionBegin; 1456 if (aij->nonew != 1) { 1457 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1458 } 1459 if (!aij->saved_values) { 1460 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1461 } 1462 1463 /* copy values over */ 1464 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1465 PetscFunctionReturn(0); 1466 } 1467 EXTERN_C_END 1468 1469 #undef __FUNC__ 1470 #define __FUNC__ /*<a name=""></a>*/"MatCreateSeqBAIJ" 1471 /*@C 1472 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1473 compressed row) format. For good matrix assembly performance the 1474 user should preallocate the matrix storage by setting the parameter nz 1475 (or the array nnz). By setting these parameters accurately, performance 1476 during matrix assembly can be increased by more than a factor of 50. 1477 1478 Collective on MPI_Comm 1479 1480 Input Parameters: 1481 + comm - MPI communicator, set to PETSC_COMM_SELF 1482 . bs - size of block 1483 . m - number of rows 1484 . n - number of columns 1485 . nz - number of block nonzeros per block row (same for all rows) 1486 - nnz - array containing the number of block nonzeros in the various block rows 1487 (possibly different for each block row) or PETSC_NULL 1488 1489 Output Parameter: 1490 . A - the matrix 1491 1492 Options Database Keys: 1493 . -mat_no_unroll - uses code that does not unroll the loops in the 1494 block calculations (much slower) 1495 . -mat_block_size - size of the blocks to use 1496 1497 Level: intermediate 1498 1499 Notes: 1500 The block AIJ format is fully compatible with standard Fortran 77 1501 storage. That is, the stored row and column indices can begin at 1502 either one (as in Fortran) or zero. See the users' manual for details. 1503 1504 Specify the preallocated storage with either nz or nnz (not both). 1505 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1506 allocation. For additional details, see the users manual chapter on 1507 matrices. 1508 1509 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1510 @*/ 1511 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1512 { 1513 Mat B; 1514 Mat_SeqBAIJ *b; 1515 int i,len,ierr,mbs,nbs,bs2,size; 1516 PetscTruth flg; 1517 1518 PetscFunctionBegin; 1519 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1520 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1521 1522 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1523 mbs = m/bs; 1524 nbs = n/bs; 1525 bs2 = bs*bs; 1526 1527 if (mbs*bs!=m || nbs*bs!=n) { 1528 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1529 } 1530 1531 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"nz cannot be less than -2: value %d",nz); 1532 if (nnz) { 1533 for (i=0; i<mbs; i++) { 1534 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1535 } 1536 } 1537 1538 *A = 0; 1539 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,"Mat",comm,MatDestroy,MatView); 1540 PLogObjectCreate(B); 1541 B->data = (void*)(b = PetscNew(Mat_SeqBAIJ));CHKPTRQ(b); 1542 ierr = PetscMemzero(b,sizeof(Mat_SeqBAIJ));CHKERRQ(ierr); 1543 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1544 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1545 if (!flg) { 1546 switch (bs) { 1547 case 1: 1548 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1549 B->ops->solve = MatSolve_SeqBAIJ_1; 1550 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 1551 B->ops->mult = MatMult_SeqBAIJ_1; 1552 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1553 break; 1554 case 2: 1555 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1556 B->ops->solve = MatSolve_SeqBAIJ_2; 1557 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 1558 B->ops->mult = MatMult_SeqBAIJ_2; 1559 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1560 break; 1561 case 3: 1562 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1563 B->ops->solve = MatSolve_SeqBAIJ_3; 1564 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3; 1565 B->ops->mult = MatMult_SeqBAIJ_3; 1566 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1567 break; 1568 case 4: 1569 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1570 B->ops->solve = MatSolve_SeqBAIJ_4; 1571 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4; 1572 B->ops->mult = MatMult_SeqBAIJ_4; 1573 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1574 break; 1575 case 5: 1576 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1577 B->ops->solve = MatSolve_SeqBAIJ_5; 1578 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5; 1579 B->ops->mult = MatMult_SeqBAIJ_5; 1580 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1581 break; 1582 case 6: 1583 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 1584 B->ops->solve = MatSolve_SeqBAIJ_6; 1585 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6; 1586 B->ops->mult = MatMult_SeqBAIJ_6; 1587 B->ops->multadd = MatMultAdd_SeqBAIJ_6; 1588 break; 1589 case 7: 1590 B->ops->mult = MatMult_SeqBAIJ_7; 1591 B->ops->solve = MatSolve_SeqBAIJ_7; 1592 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7; 1593 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1594 break; 1595 } 1596 } 1597 B->ops->destroy = MatDestroy_SeqBAIJ; 1598 B->ops->view = MatView_SeqBAIJ; 1599 B->factor = 0; 1600 B->lupivotthreshold = 1.0; 1601 B->mapping = 0; 1602 b->row = 0; 1603 b->col = 0; 1604 b->icol = 0; 1605 b->reallocs = 0; 1606 b->saved_values = 0; 1607 1608 b->m = m; B->m = m; B->M = m; 1609 b->n = n; B->n = n; B->N = n; 1610 1611 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1612 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1613 1614 b->mbs = mbs; 1615 b->nbs = nbs; 1616 b->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(b->imax); 1617 if (!nnz) { 1618 if (nz == PETSC_DEFAULT) nz = 5; 1619 else if (nz <= 0) nz = 1; 1620 for (i=0; i<mbs; i++) b->imax[i] = nz; 1621 nz = nz*mbs; 1622 } else { 1623 nz = 0; 1624 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1625 } 1626 1627 /* allocate the matrix space */ 1628 len = nz*sizeof(int) + nz*bs2*sizeof(MatScalar) + (b->m+1)*sizeof(int); 1629 b->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(b->a); 1630 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1631 b->j = (int*)(b->a + nz*bs2); 1632 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 1633 b->i = b->j + nz; 1634 b->singlemalloc = PETSC_TRUE; 1635 1636 b->i[0] = 0; 1637 for (i=1; i<mbs+1; i++) { 1638 b->i[i] = b->i[i-1] + b->imax[i-1]; 1639 } 1640 1641 /* b->ilen will count nonzeros in each block row so far. */ 1642 b->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int)); 1643 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1644 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1645 1646 b->bs = bs; 1647 b->bs2 = bs2; 1648 b->mbs = mbs; 1649 b->nz = 0; 1650 b->maxnz = nz*bs2; 1651 b->sorted = PETSC_FALSE; 1652 b->roworiented = PETSC_TRUE; 1653 b->nonew = 0; 1654 b->diag = 0; 1655 b->solve_work = 0; 1656 b->mult_work = 0; 1657 b->spptr = 0; 1658 B->info.nz_unneeded = (PetscReal)b->maxnz; 1659 b->keepzeroedrows = PETSC_FALSE; 1660 1661 *A = B; 1662 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 1663 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr); } 1664 1665 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1666 "MatStoreValues_SeqBAIJ", 1667 (void*)MatStoreValues_SeqBAIJ);CHKERRQ(ierr); 1668 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1669 "MatRetrieveValues_SeqBAIJ", 1670 (void*)MatRetrieveValues_SeqBAIJ);CHKERRQ(ierr); 1671 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1672 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1673 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1674 PetscFunctionReturn(0); 1675 } 1676 1677 #undef __FUNC__ 1678 #define __FUNC__ /*<a name=""></a>*/"MatDuplicate_SeqBAIJ" 1679 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1680 { 1681 Mat C; 1682 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ*)A->data; 1683 int i,len,mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1684 1685 PetscFunctionBegin; 1686 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1687 1688 *B = 0; 1689 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,"Mat",A->comm,MatDestroy,MatView); 1690 PLogObjectCreate(C); 1691 C->data = (void*)(c = PetscNew(Mat_SeqBAIJ));CHKPTRQ(c); 1692 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1693 C->ops->destroy = MatDestroy_SeqBAIJ; 1694 C->ops->view = MatView_SeqBAIJ; 1695 C->factor = A->factor; 1696 c->row = 0; 1697 c->col = 0; 1698 c->icol = 0; 1699 c->saved_values = 0; 1700 c->keepzeroedrows = a->keepzeroedrows; 1701 C->assembled = PETSC_TRUE; 1702 1703 c->m = C->m = a->m; 1704 c->n = C->n = a->n; 1705 C->M = a->m; 1706 C->N = a->n; 1707 1708 c->bs = a->bs; 1709 c->bs2 = a->bs2; 1710 c->mbs = a->mbs; 1711 c->nbs = a->nbs; 1712 1713 c->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->imax); 1714 c->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->ilen); 1715 for (i=0; i<mbs; i++) { 1716 c->imax[i] = a->imax[i]; 1717 c->ilen[i] = a->ilen[i]; 1718 } 1719 1720 /* allocate the matrix space */ 1721 c->singlemalloc = PETSC_TRUE; 1722 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1723 c->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(c->a); 1724 c->j = (int*)(c->a + nz*bs2); 1725 c->i = c->j + nz; 1726 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1727 if (mbs > 0) { 1728 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1729 if (cpvalues == MAT_COPY_VALUES) { 1730 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1731 } else { 1732 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1733 } 1734 } 1735 1736 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1737 c->sorted = a->sorted; 1738 c->roworiented = a->roworiented; 1739 c->nonew = a->nonew; 1740 1741 if (a->diag) { 1742 c->diag = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->diag); 1743 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1744 for (i=0; i<mbs; i++) { 1745 c->diag[i] = a->diag[i]; 1746 } 1747 } else c->diag = 0; 1748 c->nz = a->nz; 1749 c->maxnz = a->maxnz; 1750 c->solve_work = 0; 1751 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1752 c->mult_work = 0; 1753 *B = C; 1754 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1755 PetscFunctionReturn(0); 1756 } 1757 1758 #undef __FUNC__ 1759 #define __FUNC__ /*<a name=""></a>*/"MatLoad_SeqBAIJ" 1760 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1761 { 1762 Mat_SeqBAIJ *a; 1763 Mat B; 1764 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1765 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1766 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1767 int *masked,nmask,tmp,bs2,ishift; 1768 Scalar *aa; 1769 MPI_Comm comm = ((PetscObject)viewer)->comm; 1770 1771 PetscFunctionBegin; 1772 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1773 bs2 = bs*bs; 1774 1775 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1776 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1777 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1778 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1779 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1780 M = header[1]; N = header[2]; nz = header[3]; 1781 1782 if (header[3] < 0) { 1783 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1784 } 1785 1786 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1787 1788 /* 1789 This code adds extra rows to make sure the number of rows is 1790 divisible by the blocksize 1791 */ 1792 mbs = M/bs; 1793 extra_rows = bs - M + bs*(mbs); 1794 if (extra_rows == bs) extra_rows = 0; 1795 else mbs++; 1796 if (extra_rows) { 1797 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1798 } 1799 1800 /* read in row lengths */ 1801 rowlengths = (int*)PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1802 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1803 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1804 1805 /* read in column indices */ 1806 jj = (int*)PetscMalloc((nz+extra_rows)*sizeof(int));CHKPTRQ(jj); 1807 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1808 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1809 1810 /* loop over row lengths determining block row lengths */ 1811 browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1812 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1813 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1814 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1815 masked = mask + mbs; 1816 rowcount = 0; nzcount = 0; 1817 for (i=0; i<mbs; i++) { 1818 nmask = 0; 1819 for (j=0; j<bs; j++) { 1820 kmax = rowlengths[rowcount]; 1821 for (k=0; k<kmax; k++) { 1822 tmp = jj[nzcount++]/bs; 1823 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1824 } 1825 rowcount++; 1826 } 1827 browlengths[i] += nmask; 1828 /* zero out the mask elements we set */ 1829 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1830 } 1831 1832 /* create our matrix */ 1833 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1834 B = *A; 1835 a = (Mat_SeqBAIJ*)B->data; 1836 1837 /* set matrix "i" values */ 1838 a->i[0] = 0; 1839 for (i=1; i<= mbs; i++) { 1840 a->i[i] = a->i[i-1] + browlengths[i-1]; 1841 a->ilen[i-1] = browlengths[i-1]; 1842 } 1843 a->nz = 0; 1844 for (i=0; i<mbs; i++) a->nz += browlengths[i]; 1845 1846 /* read in nonzero values */ 1847 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1848 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1849 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1850 1851 /* set "a" and "j" values into matrix */ 1852 nzcount = 0; jcount = 0; 1853 for (i=0; i<mbs; i++) { 1854 nzcountb = nzcount; 1855 nmask = 0; 1856 for (j=0; j<bs; j++) { 1857 kmax = rowlengths[i*bs+j]; 1858 for (k=0; k<kmax; k++) { 1859 tmp = jj[nzcount++]/bs; 1860 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1861 } 1862 rowcount++; 1863 } 1864 /* sort the masked values */ 1865 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1866 1867 /* set "j" values into matrix */ 1868 maskcount = 1; 1869 for (j=0; j<nmask; j++) { 1870 a->j[jcount++] = masked[j]; 1871 mask[masked[j]] = maskcount++; 1872 } 1873 /* set "a" values into matrix */ 1874 ishift = bs2*a->i[i]; 1875 for (j=0; j<bs; j++) { 1876 kmax = rowlengths[i*bs+j]; 1877 for (k=0; k<kmax; k++) { 1878 tmp = jj[nzcountb]/bs ; 1879 block = mask[tmp] - 1; 1880 point = jj[nzcountb] - bs*tmp; 1881 idx = ishift + bs2*block + j + bs*point; 1882 a->a[idx] = (MatScalar)aa[nzcountb++]; 1883 } 1884 } 1885 /* zero out the mask elements we set */ 1886 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1887 } 1888 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1889 1890 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1891 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1892 ierr = PetscFree(aa);CHKERRQ(ierr); 1893 ierr = PetscFree(jj);CHKERRQ(ierr); 1894 ierr = PetscFree(mask);CHKERRQ(ierr); 1895 1896 B->assembled = PETSC_TRUE; 1897 1898 ierr = MatView_Private(B);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 1903 1904