1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.127 1998/03/13 19:06:03 balay Exp balay $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the BAIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/baij/seq/baij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "petsc.h" 16 17 #define CHUNKSIZE 10 18 19 #undef __FUNC__ 20 #define __FUNC__ "MatMarkDiag_SeqBAIJ" 21 int MatMarkDiag_SeqBAIJ(Mat A) 22 { 23 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 24 int i,j, *diag, m = a->mbs; 25 26 PetscFunctionBegin; 27 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 28 PLogObjectMemory(A,(m+1)*sizeof(int)); 29 for ( i=0; i<m; i++ ) { 30 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 31 if (a->j[j] == i) { 32 diag[i] = j; 33 break; 34 } 35 } 36 } 37 a->diag = diag; 38 PetscFunctionReturn(0); 39 } 40 41 42 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 43 44 #undef __FUNC__ 45 #define __FUNC__ "MatGetRowIJ_SeqBAIJ" 46 static int MatGetRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 47 PetscTruth *done) 48 { 49 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 50 int ierr, n = a->mbs,i; 51 52 PetscFunctionBegin; 53 *nn = n; 54 if (!ia) PetscFunctionReturn(0); 55 if (symmetric) { 56 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 57 } else if (oshift == 1) { 58 /* temporarily add 1 to i and j indices */ 59 int nz = a->i[n] + 1; 60 for ( i=0; i<nz; i++ ) a->j[i]++; 61 for ( i=0; i<n+1; i++ ) a->i[i]++; 62 *ia = a->i; *ja = a->j; 63 } else { 64 *ia = a->i; *ja = a->j; 65 } 66 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNC__ 71 #define __FUNC__ "MatRestoreRowIJ_SeqBAIJ" 72 static int MatRestoreRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 73 PetscTruth *done) 74 { 75 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 76 int i,n = a->mbs; 77 78 PetscFunctionBegin; 79 if (!ia) PetscFunctionReturn(0); 80 if (symmetric) { 81 PetscFree(*ia); 82 PetscFree(*ja); 83 } else if (oshift == 1) { 84 int nz = a->i[n]; 85 for ( i=0; i<nz; i++ ) a->j[i]--; 86 for ( i=0; i<n+1; i++ ) a->i[i]--; 87 } 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNC__ 92 #define __FUNC__ "MatGetBlockSize_SeqBAIJ" 93 int MatGetBlockSize_SeqBAIJ(Mat mat, int *bs) 94 { 95 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *) mat->data; 96 97 PetscFunctionBegin; 98 *bs = baij->bs; 99 PetscFunctionReturn(0); 100 } 101 102 103 #undef __FUNC__ 104 #define __FUNC__ "MatDestroy_SeqBAIJ" 105 int MatDestroy_SeqBAIJ(PetscObject obj) 106 { 107 Mat A = (Mat) obj; 108 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 109 int ierr; 110 111 #if defined(USE_PETSC_LOG) 112 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 113 #endif 114 PetscFree(a->a); 115 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 116 if (a->diag) PetscFree(a->diag); 117 if (a->ilen) PetscFree(a->ilen); 118 if (a->imax) PetscFree(a->imax); 119 if (a->solve_work) PetscFree(a->solve_work); 120 if (a->mult_work) PetscFree(a->mult_work); 121 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 122 PetscFree(a); 123 PLogObjectDestroy(A); 124 PetscHeaderDestroy(A); 125 PetscFunctionReturn(0); 126 } 127 128 #undef __FUNC__ 129 #define __FUNC__ "MatSetOption_SeqBAIJ" 130 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 131 { 132 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 133 134 PetscFunctionBegin; 135 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 136 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 137 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 138 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 139 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 140 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 141 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 142 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 143 else if (op == MAT_ROWS_SORTED || 144 op == MAT_ROWS_UNSORTED || 145 op == MAT_SYMMETRIC || 146 op == MAT_STRUCTURALLY_SYMMETRIC || 147 op == MAT_YES_NEW_DIAGONALS || 148 op == MAT_IGNORE_OFF_PROC_ENTRIES || 149 op == MAT_USE_HASH_TABLE) { 150 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 151 } else if (op == MAT_NO_NEW_DIAGONALS) { 152 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 153 } else { 154 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 155 } 156 PetscFunctionReturn(0); 157 } 158 159 160 #undef __FUNC__ 161 #define __FUNC__ "MatGetSize_SeqBAIJ" 162 int MatGetSize_SeqBAIJ(Mat A,int *m,int *n) 163 { 164 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 165 166 PetscFunctionBegin; 167 if (m) *m = a->m; 168 if (n) *n = a->n; 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNC__ 173 #define __FUNC__ "MatGetOwnershipRange_SeqBAIJ" 174 int MatGetOwnershipRange_SeqBAIJ(Mat A,int *m,int *n) 175 { 176 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 177 178 PetscFunctionBegin; 179 *m = 0; *n = a->m; 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNC__ 184 #define __FUNC__ "MatGetRow_SeqBAIJ" 185 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 186 { 187 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 188 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2; 189 Scalar *aa,*v_i,*aa_i; 190 191 PetscFunctionBegin; 192 bs = a->bs; 193 ai = a->i; 194 aj = a->j; 195 aa = a->a; 196 bs2 = a->bs2; 197 198 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 199 200 bn = row/bs; /* Block number */ 201 bp = row % bs; /* Block Position */ 202 M = ai[bn+1] - ai[bn]; 203 *nz = bs*M; 204 205 if (v) { 206 *v = 0; 207 if (*nz) { 208 *v = (Scalar *) PetscMalloc( (*nz)*sizeof(Scalar) ); CHKPTRQ(*v); 209 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 210 v_i = *v + i*bs; 211 aa_i = aa + bs2*(ai[bn] + i); 212 for ( j=bp,k=0; j<bs2; j+=bs,k++ ) {v_i[k] = aa_i[j];} 213 } 214 } 215 } 216 217 if (idx) { 218 *idx = 0; 219 if (*nz) { 220 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 221 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 222 idx_i = *idx + i*bs; 223 itmp = bs*aj[ai[bn] + i]; 224 for ( j=0; j<bs; j++ ) {idx_i[j] = itmp++;} 225 } 226 } 227 } 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNC__ 232 #define __FUNC__ "MatRestoreRow_SeqBAIJ" 233 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 234 { 235 PetscFunctionBegin; 236 if (idx) {if (*idx) PetscFree(*idx);} 237 if (v) {if (*v) PetscFree(*v);} 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNC__ 242 #define __FUNC__ "MatTranspose_SeqBAIJ" 243 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 244 { 245 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 246 Mat C; 247 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 248 int *rows,*cols,bs2=a->bs2; 249 Scalar *array=a->a; 250 251 PetscFunctionBegin; 252 if (B==PETSC_NULL && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 253 col = (int *) PetscMalloc((1+nbs)*sizeof(int)); CHKPTRQ(col); 254 PetscMemzero(col,(1+nbs)*sizeof(int)); 255 256 for ( i=0; i<ai[mbs]; i++ ) col[aj[i]] += 1; 257 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C); CHKERRQ(ierr); 258 PetscFree(col); 259 rows = (int *) PetscMalloc(2*bs*sizeof(int)); CHKPTRQ(rows); 260 cols = rows + bs; 261 for ( i=0; i<mbs; i++ ) { 262 cols[0] = i*bs; 263 for (k=1; k<bs; k++ ) cols[k] = cols[k-1] + 1; 264 len = ai[i+1] - ai[i]; 265 for ( j=0; j<len; j++ ) { 266 rows[0] = (*aj++)*bs; 267 for (k=1; k<bs; k++ ) rows[k] = rows[k-1] + 1; 268 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES); CHKERRQ(ierr); 269 array += bs2; 270 } 271 } 272 PetscFree(rows); 273 274 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 275 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 276 277 if (B != PETSC_NULL) { 278 *B = C; 279 } else { 280 PetscOps *Abops; 281 struct _MatOps *Aops; 282 283 /* This isn't really an in-place transpose */ 284 PetscFree(a->a); 285 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 286 if (a->diag) PetscFree(a->diag); 287 if (a->ilen) PetscFree(a->ilen); 288 if (a->imax) PetscFree(a->imax); 289 if (a->solve_work) PetscFree(a->solve_work); 290 PetscFree(a); 291 292 /* 293 This is horrible, horrible code. We need to keep the 294 A pointers for the bops and ops but copy everything 295 else from C. 296 */ 297 Abops = A->bops; 298 Aops = A->ops; 299 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 300 A->bops = Abops; 301 A->ops = Aops; 302 303 PetscHeaderDestroy(C); 304 } 305 PetscFunctionReturn(0); 306 } 307 308 309 310 311 #undef __FUNC__ 312 #define __FUNC__ "MatView_SeqBAIJ_Binary" 313 static int MatView_SeqBAIJ_Binary(Mat A,Viewer viewer) 314 { 315 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 316 int i, fd, *col_lens, ierr, bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 317 Scalar *aa; 318 319 PetscFunctionBegin; 320 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 321 col_lens = (int *) PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 322 col_lens[0] = MAT_COOKIE; 323 324 col_lens[1] = a->m; 325 col_lens[2] = a->n; 326 col_lens[3] = a->nz*bs2; 327 328 /* store lengths of each row and write (including header) to file */ 329 count = 0; 330 for ( i=0; i<a->mbs; i++ ) { 331 for ( j=0; j<bs; j++ ) { 332 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 333 } 334 } 335 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 336 PetscFree(col_lens); 337 338 /* store column indices (zero start index) */ 339 jj = (int *) PetscMalloc( (a->nz+1)*bs2*sizeof(int) ); CHKPTRQ(jj); 340 count = 0; 341 for ( i=0; i<a->mbs; i++ ) { 342 for ( j=0; j<bs; j++ ) { 343 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 344 for ( l=0; l<bs; l++ ) { 345 jj[count++] = bs*a->j[k] + l; 346 } 347 } 348 } 349 } 350 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0); CHKERRQ(ierr); 351 PetscFree(jj); 352 353 /* store nonzero values */ 354 aa = (Scalar *) PetscMalloc((a->nz+1)*bs2*sizeof(Scalar)); CHKPTRQ(aa); 355 count = 0; 356 for ( i=0; i<a->mbs; i++ ) { 357 for ( j=0; j<bs; j++ ) { 358 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 359 for ( l=0; l<bs; l++ ) { 360 aa[count++] = a->a[bs2*k + l*bs + j]; 361 } 362 } 363 } 364 } 365 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 366 PetscFree(aa); 367 PetscFunctionReturn(0); 368 } 369 370 #undef __FUNC__ 371 #define __FUNC__ "MatView_SeqBAIJ_ASCII" 372 static int MatView_SeqBAIJ_ASCII(Mat A,Viewer viewer) 373 { 374 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 375 int ierr, i,j,format,bs = a->bs,k,l,bs2=a->bs2; 376 FILE *fd; 377 char *outputname; 378 379 PetscFunctionBegin; 380 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 381 ierr = ViewerFileGetOutputname_Private(viewer,&outputname);CHKERRQ(ierr); 382 ierr = ViewerGetFormat(viewer,&format); 383 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 384 fprintf(fd," block size is %d\n",bs); 385 } 386 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 387 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 388 } 389 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 390 for ( i=0; i<a->mbs; i++ ) { 391 for ( j=0; j<bs; j++ ) { 392 fprintf(fd,"row %d:",i*bs+j); 393 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 394 for ( l=0; l<bs; l++ ) { 395 #if defined(USE_PETSC_COMPLEX) 396 if (imag(a->a[bs2*k + l*bs + j]) > 0.0 && real(a->a[bs2*k + l*bs + j]) != 0.0) 397 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 398 real(a->a[bs2*k + l*bs + j]),imag(a->a[bs2*k + l*bs + j])); 399 else if (imag(a->a[bs2*k + l*bs + j]) < 0.0 && real(a->a[bs2*k + l*bs + j]) != 0.0) 400 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 401 real(a->a[bs2*k + l*bs + j]),-imag(a->a[bs2*k + l*bs + j])); 402 else if (real(a->a[bs2*k + l*bs + j]) != 0.0) 403 fprintf(fd," %d %g ",bs*a->j[k]+l,real(a->a[bs2*k + l*bs + j])); 404 #else 405 if (a->a[bs2*k + l*bs + j] != 0.0) 406 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 407 #endif 408 } 409 } 410 fprintf(fd,"\n"); 411 } 412 } 413 } 414 else { 415 for ( i=0; i<a->mbs; i++ ) { 416 for ( j=0; j<bs; j++ ) { 417 fprintf(fd,"row %d:",i*bs+j); 418 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 419 for ( l=0; l<bs; l++ ) { 420 #if defined(USE_PETSC_COMPLEX) 421 if (imag(a->a[bs2*k + l*bs + j]) > 0.0) { 422 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 423 real(a->a[bs2*k + l*bs + j]),imag(a->a[bs2*k + l*bs + j])); 424 } 425 else if (imag(a->a[bs2*k + l*bs + j]) < 0.0) { 426 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 427 real(a->a[bs2*k + l*bs + j]),-imag(a->a[bs2*k + l*bs + j])); 428 } 429 else { 430 fprintf(fd," %d %g ",bs*a->j[k]+l,real(a->a[bs2*k + l*bs + j])); 431 } 432 #else 433 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 434 #endif 435 } 436 } 437 fprintf(fd,"\n"); 438 } 439 } 440 } 441 fflush(fd); 442 PetscFunctionReturn(0); 443 } 444 445 #undef __FUNC__ 446 #define __FUNC__ "MatView_SeqBAIJ_Draw" 447 static int MatView_SeqBAIJ_Draw(Mat A,Viewer viewer) 448 { 449 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 450 int row,ierr,i,j,k,l,mbs=a->mbs,pause,color,bs=a->bs,bs2=a->bs2; 451 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r; 452 Scalar *aa; 453 Draw draw; 454 DrawButton button; 455 PetscTruth isnull; 456 457 PetscFunctionBegin; 458 ierr = ViewerDrawGetDraw(viewer,&draw);CHKERRQ(ierr); 459 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 460 461 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 462 xr += w; yr += h; xl = -w; yl = -h; 463 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 464 /* loop over matrix elements drawing boxes */ 465 color = DRAW_BLUE; 466 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 467 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 468 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 469 x_l = a->j[j]*bs; x_r = x_l + 1.0; 470 aa = a->a + j*bs2; 471 for ( k=0; k<bs; k++ ) { 472 for ( l=0; l<bs; l++ ) { 473 if (PetscReal(*aa++) >= 0.) continue; 474 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 475 } 476 } 477 } 478 } 479 color = DRAW_CYAN; 480 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 481 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 482 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 483 x_l = a->j[j]*bs; x_r = x_l + 1.0; 484 aa = a->a + j*bs2; 485 for ( k=0; k<bs; k++ ) { 486 for ( l=0; l<bs; l++ ) { 487 if (PetscReal(*aa++) != 0.) continue; 488 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 489 } 490 } 491 } 492 } 493 494 color = DRAW_RED; 495 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 496 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 497 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 498 x_l = a->j[j]*bs; x_r = x_l + 1.0; 499 aa = a->a + j*bs2; 500 for ( k=0; k<bs; k++ ) { 501 for ( l=0; l<bs; l++ ) { 502 if (PetscReal(*aa++) <= 0.) continue; 503 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 504 } 505 } 506 } 507 } 508 509 DrawSynchronizedFlush(draw); 510 DrawGetPause(draw,&pause); 511 if (pause >= 0) { PetscSleep(pause); PetscFunctionReturn(0);} 512 513 /* allow the matrix to zoom or shrink */ 514 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); 515 while (button != BUTTON_RIGHT) { 516 DrawSynchronizedClear(draw); 517 if (button == BUTTON_LEFT) scale = .5; 518 else if (button == BUTTON_CENTER) scale = 2.; 519 xl = scale*(xl + w - xc) + xc - w*scale; 520 xr = scale*(xr - w - xc) + xc + w*scale; 521 yl = scale*(yl + h - yc) + yc - h*scale; 522 yr = scale*(yr - h - yc) + yc + h*scale; 523 w *= scale; h *= scale; 524 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 525 color = DRAW_BLUE; 526 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 527 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 528 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 529 x_l = a->j[j]*bs; x_r = x_l + 1.0; 530 aa = a->a + j*bs2; 531 for ( k=0; k<bs; k++ ) { 532 for ( l=0; l<bs; l++ ) { 533 if (PetscReal(*aa++) >= 0.) continue; 534 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 535 } 536 } 537 } 538 } 539 color = DRAW_CYAN; 540 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 541 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 542 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 543 x_l = a->j[j]*bs; x_r = x_l + 1.0; 544 aa = a->a + j*bs2; 545 for ( k=0; k<bs; k++ ) { 546 for ( l=0; l<bs; l++ ) { 547 if (PetscReal(*aa++) != 0.) continue; 548 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 549 } 550 } 551 } 552 } 553 554 color = DRAW_RED; 555 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 556 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 557 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 558 x_l = a->j[j]*bs; x_r = x_l + 1.0; 559 aa = a->a + j*bs2; 560 for ( k=0; k<bs; k++ ) { 561 for ( l=0; l<bs; l++ ) { 562 if (PetscReal(*aa++) <= 0.) continue; 563 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 564 } 565 } 566 } 567 } 568 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); 569 } 570 PetscFunctionReturn(0); 571 } 572 573 #undef __FUNC__ 574 #define __FUNC__ "MatView_SeqBAIJ" 575 int MatView_SeqBAIJ(PetscObject obj,Viewer viewer) 576 { 577 Mat A = (Mat) obj; 578 ViewerType vtype; 579 int ierr; 580 581 PetscFunctionBegin; 582 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 583 if (vtype == MATLAB_VIEWER) { 584 SETERRQ(PETSC_ERR_SUP,0,"Matlab viewer not supported"); 585 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 586 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 587 } else if (vtype == BINARY_FILE_VIEWER) { 588 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 589 } else if (vtype == DRAW_VIEWER) { 590 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 591 } 592 PetscFunctionReturn(0); 593 } 594 595 596 #undef __FUNC__ 597 #define __FUNC__ "MatGetValues_SeqBAIJ" 598 int MatGetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 599 { 600 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 601 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 602 int *ai = a->i, *ailen = a->ilen; 603 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 604 Scalar *ap, *aa = a->a, zero = 0.0; 605 606 PetscFunctionBegin; 607 for ( k=0; k<m; k++ ) { /* loop over rows */ 608 row = im[k]; brow = row/bs; 609 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 610 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 611 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 612 nrow = ailen[brow]; 613 for ( l=0; l<n; l++ ) { /* loop over columns */ 614 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 615 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 616 col = in[l] ; 617 bcol = col/bs; 618 cidx = col%bs; 619 ridx = row%bs; 620 high = nrow; 621 low = 0; /* assume unsorted */ 622 while (high-low > 5) { 623 t = (low+high)/2; 624 if (rp[t] > bcol) high = t; 625 else low = t; 626 } 627 for ( i=low; i<high; i++ ) { 628 if (rp[i] > bcol) break; 629 if (rp[i] == bcol) { 630 *v++ = ap[bs2*i+bs*cidx+ridx]; 631 goto finished; 632 } 633 } 634 *v++ = zero; 635 finished:; 636 } 637 } 638 PetscFunctionReturn(0); 639 } 640 641 642 #undef __FUNC__ 643 #define __FUNC__ "MatSetValuesBlocked_SeqBAIJ" 644 int MatSetValuesBlocked_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 645 { 646 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 647 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 648 int *imax=a->imax,*ai=a->i,*ailen=a->ilen, roworiented=a->roworiented; 649 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval; 650 Scalar *ap,*value=v,*aa=a->a,*bap; 651 652 PetscFunctionBegin; 653 if (roworiented) { 654 stepval = (n-1)*bs; 655 } else { 656 stepval = (m-1)*bs; 657 } 658 for ( k=0; k<m; k++ ) { /* loop over added rows */ 659 row = im[k]; 660 #if defined(USE_PETSC_BOPT_g) 661 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 662 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 663 #endif 664 rp = aj + ai[row]; 665 ap = aa + bs2*ai[row]; 666 rmax = imax[row]; 667 nrow = ailen[row]; 668 low = 0; 669 for ( l=0; l<n; l++ ) { /* loop over added columns */ 670 #if defined(USE_PETSC_BOPT_g) 671 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 672 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 673 #endif 674 col = in[l]; 675 if (roworiented) { 676 value = v + k*(stepval+bs)*bs + l*bs; 677 } else { 678 value = v + l*(stepval+bs)*bs + k*bs; 679 } 680 if (!sorted) low = 0; high = nrow; 681 while (high-low > 7) { 682 t = (low+high)/2; 683 if (rp[t] > col) high = t; 684 else low = t; 685 } 686 for ( i=low; i<high; i++ ) { 687 if (rp[i] > col) break; 688 if (rp[i] == col) { 689 bap = ap + bs2*i; 690 if (roworiented) { 691 if (is == ADD_VALUES) { 692 for ( ii=0; ii<bs; ii++,value+=stepval ) { 693 for (jj=ii; jj<bs2; jj+=bs ) { 694 bap[jj] += *value++; 695 } 696 } 697 } else { 698 for ( ii=0; ii<bs; ii++,value+=stepval ) { 699 for (jj=ii; jj<bs2; jj+=bs ) { 700 bap[jj] = *value++; 701 } 702 } 703 } 704 } else { 705 if (is == ADD_VALUES) { 706 for ( ii=0; ii<bs; ii++,value+=stepval ) { 707 for (jj=0; jj<bs; jj++ ) { 708 *bap++ += *value++; 709 } 710 } 711 } else { 712 for ( ii=0; ii<bs; ii++,value+=stepval ) { 713 for (jj=0; jj<bs; jj++ ) { 714 *bap++ = *value++; 715 } 716 } 717 } 718 } 719 goto noinsert2; 720 } 721 } 722 if (nonew == 1) goto noinsert2; 723 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 724 if (nrow >= rmax) { 725 /* there is no extra room in row, therefore enlarge */ 726 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 727 Scalar *new_a; 728 729 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 730 731 /* malloc new storage space */ 732 len = new_nz*(sizeof(int)+bs2*sizeof(Scalar))+(a->mbs+1)*sizeof(int); 733 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 734 new_j = (int *) (new_a + bs2*new_nz); 735 new_i = new_j + new_nz; 736 737 /* copy over old data into new slots */ 738 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 739 for ( ii=row+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 740 PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int)); 741 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 742 PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int)); 743 PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(Scalar)); 744 PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(Scalar)); 745 PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(Scalar)); 746 /* free up old matrix storage */ 747 PetscFree(a->a); 748 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 749 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 750 a->singlemalloc = 1; 751 752 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 753 rmax = imax[row] = imax[row] + CHUNKSIZE; 754 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(Scalar))); 755 a->maxnz += bs2*CHUNKSIZE; 756 a->reallocs++; 757 a->nz++; 758 } 759 N = nrow++ - 1; 760 /* shift up all the later entries in this row */ 761 for ( ii=N; ii>=i; ii-- ) { 762 rp[ii+1] = rp[ii]; 763 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(Scalar)); 764 } 765 if (N >= i) PetscMemzero(ap+bs2*i,bs2*sizeof(Scalar)); 766 rp[i] = col; 767 bap = ap + bs2*i; 768 if (roworiented) { 769 for ( ii=0; ii<bs; ii++,value+=stepval) { 770 for (jj=ii; jj<bs2; jj+=bs ) { 771 bap[jj] = *value++; 772 } 773 } 774 } else { 775 for ( ii=0; ii<bs; ii++,value+=stepval ) { 776 for (jj=0; jj<bs; jj++ ) { 777 *bap++ = *value++; 778 } 779 } 780 } 781 noinsert2:; 782 low = i; 783 } 784 ailen[row] = nrow; 785 } 786 PetscFunctionReturn(0); 787 } 788 789 790 #undef __FUNC__ 791 #define __FUNC__ "MatAssemblyEnd_SeqBAIJ" 792 int MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode) 793 { 794 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 795 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax; 796 int m = a->m,*ip, N, *ailen = a->ilen; 797 int mbs = a->mbs, bs2 = a->bs2,rmax = 0; 798 Scalar *aa = a->a, *ap; 799 800 PetscFunctionBegin; 801 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 802 803 if (m) rmax = ailen[0]; 804 for ( i=1; i<mbs; i++ ) { 805 /* move each row back by the amount of empty slots (fshift) before it*/ 806 fshift += imax[i-1] - ailen[i-1]; 807 rmax = PetscMax(rmax,ailen[i]); 808 if (fshift) { 809 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 810 N = ailen[i]; 811 for ( j=0; j<N; j++ ) { 812 ip[j-fshift] = ip[j]; 813 PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(Scalar)); 814 } 815 } 816 ai[i] = ai[i-1] + ailen[i-1]; 817 } 818 if (mbs) { 819 fshift += imax[mbs-1] - ailen[mbs-1]; 820 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 821 } 822 /* reset ilen and imax for each row */ 823 for ( i=0; i<mbs; i++ ) { 824 ailen[i] = imax[i] = ai[i+1] - ai[i]; 825 } 826 a->nz = ai[mbs]; 827 828 /* diagonals may have moved, so kill the diagonal pointers */ 829 if (fshift && a->diag) { 830 PetscFree(a->diag); 831 PLogObjectMemory(A,-(m+1)*sizeof(int)); 832 a->diag = 0; 833 } 834 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 835 m,a->n,a->bs,fshift*bs2,a->nz*bs2); 836 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Number of mallocs during MatSetValues is %d\n", 837 a->reallocs); 838 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 839 a->reallocs = 0; 840 A->info.nz_unneeded = (double)fshift*bs2; 841 842 PetscFunctionReturn(0); 843 } 844 845 846 /* idx should be of length atlease bs */ 847 #undef __FUNC__ 848 #define __FUNC__ "MatZeroRows_SeqBAIJ_Check_Block" 849 static int MatZeroRows_SeqBAIJ_Check_Block(int *idx, int bs, PetscTruth *flg) 850 { 851 int i,row; 852 853 PetscFunctionBegin; 854 row = idx[0]; 855 if (row%bs!=0) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 856 857 for ( i=1; i<bs; i++ ) { 858 if (row+i != idx[i]) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 859 } 860 *flg = PETSC_TRUE; 861 PetscFunctionReturn(0); 862 } 863 864 #undef __FUNC__ 865 #define __FUNC__ "MatZeroRows_SeqBAIJ" 866 int MatZeroRows_SeqBAIJ(Mat A,IS is, Scalar *diag) 867 { 868 Mat_SeqBAIJ *baij=(Mat_SeqBAIJ*)A->data; 869 IS is_local; 870 int ierr,i,j,count,m=baij->m,is_n,*is_idx,*rows,bs=baij->bs,bs2=baij->bs2; 871 PetscTruth flg; 872 Scalar zero = 0.0,*aa; 873 874 PetscFunctionBegin; 875 /* Make a copy of the IS and sort it */ 876 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 877 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 878 ierr = ISCreateGeneral(A->comm,is_n,is_idx,&is_local); CHKERRQ(ierr); 879 ierr = ISSort(is_local); CHKERRQ(ierr); 880 ierr = ISGetIndices(is_local,&rows); CHKERRQ(ierr); 881 882 i = 0; 883 while (i < is_n) { 884 if (rows[i]<0 || rows[i]>m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 885 flg = PETSC_FALSE; 886 if (i+bs <= is_n) {ierr = MatZeroRows_SeqBAIJ_Check_Block(rows+i,bs,&flg); CHKERRQ(ierr); } 887 count = (baij->i[rows[i]/bs +1] - baij->i[rows[i]/bs])*bs; 888 aa = baij->a + baij->i[rows[i]/bs]*bs2 + (rows[i]%bs); 889 if (flg) { /* There exists a block of rows to be Zerowed */ 890 if (baij->ilen[rows[i]/bs] > 0) { 891 PetscMemzero(aa,count*bs*sizeof(Scalar)); 892 baij->ilen[rows[i]/bs] = 1; 893 baij->j[baij->i[rows[i]/bs]] = rows[i]/bs; 894 } 895 i += bs; 896 } else { /* Zero out only the requested row */ 897 for ( j=0; j<count; j++ ) { 898 aa[0] = zero; 899 aa+=bs; 900 } 901 i++; 902 } 903 } 904 if (diag) { 905 for ( j=0; j<is_n; j++ ) { 906 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 907 /* ierr = MatSetValues(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); */ 908 } 909 } 910 ierr = ISRestoreIndices(is,&is_idx); CHKERRQ(ierr); 911 ierr = ISRestoreIndices(is_local,&rows); CHKERRQ(ierr); 912 ierr = ISDestroy(is_local); CHKERRQ(ierr); 913 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNC__ 918 #define __FUNC__ "MatSetValues_SeqBAIJ" 919 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 920 { 921 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 922 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 923 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 924 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 925 int ridx,cidx,bs2=a->bs2; 926 Scalar *ap,value,*aa=a->a,*bap; 927 928 PetscFunctionBegin; 929 for ( k=0; k<m; k++ ) { /* loop over added rows */ 930 row = im[k]; brow = row/bs; 931 #if defined(USE_PETSC_BOPT_g) 932 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 933 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 934 #endif 935 rp = aj + ai[brow]; 936 ap = aa + bs2*ai[brow]; 937 rmax = imax[brow]; 938 nrow = ailen[brow]; 939 low = 0; 940 for ( l=0; l<n; l++ ) { /* loop over added columns */ 941 #if defined(USE_PETSC_BOPT_g) 942 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 943 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 944 #endif 945 col = in[l]; bcol = col/bs; 946 ridx = row % bs; cidx = col % bs; 947 if (roworiented) { 948 value = *v++; 949 } else { 950 value = v[k + l*m]; 951 } 952 if (!sorted) low = 0; high = nrow; 953 while (high-low > 7) { 954 t = (low+high)/2; 955 if (rp[t] > bcol) high = t; 956 else low = t; 957 } 958 for ( i=low; i<high; i++ ) { 959 if (rp[i] > bcol) break; 960 if (rp[i] == bcol) { 961 bap = ap + bs2*i + bs*cidx + ridx; 962 if (is == ADD_VALUES) *bap += value; 963 else *bap = value; 964 goto noinsert1; 965 } 966 } 967 if (nonew == 1) goto noinsert1; 968 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 969 if (nrow >= rmax) { 970 /* there is no extra room in row, therefore enlarge */ 971 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 972 Scalar *new_a; 973 974 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 975 976 /* Malloc new storage space */ 977 len = new_nz*(sizeof(int)+bs2*sizeof(Scalar))+(a->mbs+1)*sizeof(int); 978 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 979 new_j = (int *) (new_a + bs2*new_nz); 980 new_i = new_j + new_nz; 981 982 /* copy over old data into new slots */ 983 for ( ii=0; ii<brow+1; ii++ ) {new_i[ii] = ai[ii];} 984 for ( ii=brow+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 985 PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int)); 986 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 987 PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow, 988 len*sizeof(int)); 989 PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(Scalar)); 990 PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(Scalar)); 991 PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE), 992 aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(Scalar)); 993 /* free up old matrix storage */ 994 PetscFree(a->a); 995 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 996 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 997 a->singlemalloc = 1; 998 999 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1000 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1001 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(Scalar))); 1002 a->maxnz += bs2*CHUNKSIZE; 1003 a->reallocs++; 1004 a->nz++; 1005 } 1006 N = nrow++ - 1; 1007 /* shift up all the later entries in this row */ 1008 for ( ii=N; ii>=i; ii-- ) { 1009 rp[ii+1] = rp[ii]; 1010 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(Scalar)); 1011 } 1012 if (N>=i) PetscMemzero(ap+bs2*i,bs2*sizeof(Scalar)); 1013 rp[i] = bcol; 1014 ap[bs2*i + bs*cidx + ridx] = value; 1015 noinsert1:; 1016 low = i; 1017 } 1018 ailen[brow] = nrow; 1019 } 1020 PetscFunctionReturn(0); 1021 } 1022 1023 extern int MatLUFactorSymbolic_SeqBAIJ(Mat,IS,IS,double,Mat*); 1024 extern int MatLUFactor_SeqBAIJ(Mat,IS,IS,double); 1025 extern int MatIncreaseOverlap_SeqBAIJ(Mat,int,IS*,int); 1026 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,int,MatGetSubMatrixCall,Mat*); 1027 extern int MatGetSubMatrices_SeqBAIJ(Mat,int,IS*,IS*,MatGetSubMatrixCall,Mat**); 1028 extern int MatMultTrans_SeqBAIJ(Mat,Vec,Vec); 1029 extern int MatMultTransAdd_SeqBAIJ(Mat,Vec,Vec,Vec); 1030 extern int MatScale_SeqBAIJ(Scalar*,Mat); 1031 extern int MatNorm_SeqBAIJ(Mat,NormType,double *); 1032 extern int MatEqual_SeqBAIJ(Mat,Mat, PetscTruth*); 1033 extern int MatGetDiagonal_SeqBAIJ(Mat,Vec); 1034 extern int MatDiagonalScale_SeqBAIJ(Mat,Vec,Vec); 1035 extern int MatGetInfo_SeqBAIJ(Mat,MatInfoType,MatInfo *); 1036 extern int MatZeroEntries_SeqBAIJ(Mat); 1037 1038 extern int MatSolve_SeqBAIJ_N(Mat,Vec,Vec); 1039 extern int MatSolve_SeqBAIJ_1(Mat,Vec,Vec); 1040 extern int MatSolve_SeqBAIJ_2(Mat,Vec,Vec); 1041 extern int MatSolve_SeqBAIJ_3(Mat,Vec,Vec); 1042 extern int MatSolve_SeqBAIJ_4(Mat,Vec,Vec); 1043 extern int MatSolve_SeqBAIJ_5(Mat,Vec,Vec); 1044 extern int MatSolve_SeqBAIJ_7(Mat,Vec,Vec); 1045 1046 extern int MatLUFactorNumeric_SeqBAIJ_N(Mat,Mat*); 1047 extern int MatLUFactorNumeric_SeqBAIJ_1(Mat,Mat*); 1048 extern int MatLUFactorNumeric_SeqBAIJ_2(Mat,Mat*); 1049 extern int MatLUFactorNumeric_SeqBAIJ_3(Mat,Mat*); 1050 extern int MatLUFactorNumeric_SeqBAIJ_4(Mat,Mat*); 1051 extern int MatLUFactorNumeric_SeqBAIJ_5(Mat,Mat*); 1052 extern int MatSolve_SeqBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 1053 1054 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1055 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1056 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1057 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1058 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1059 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1060 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1061 1062 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1063 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1064 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1065 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1066 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1067 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1068 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1069 1070 /* 1071 note: This can only work for identity for row and col. It would 1072 be good to check this and otherwise generate an error. 1073 */ 1074 #undef __FUNC__ 1075 #define __FUNC__ "MatILUFactor_SeqBAIJ" 1076 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,double efill,int fill) 1077 { 1078 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) inA->data; 1079 Mat outA; 1080 int ierr; 1081 1082 PetscFunctionBegin; 1083 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1084 1085 outA = inA; 1086 inA->factor = FACTOR_LU; 1087 a->row = row; 1088 a->col = col; 1089 1090 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1091 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1092 1093 if (!a->solve_work) { 1094 a->solve_work = (Scalar *) PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1095 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1096 } 1097 1098 if (!a->diag) { 1099 ierr = MatMarkDiag_SeqBAIJ(inA); CHKERRQ(ierr); 1100 } 1101 ierr = MatLUFactorNumeric(inA,&outA); CHKERRQ(ierr); 1102 1103 /* 1104 Blocksize 4 has a special faster solver for ILU(0) factorization 1105 with natural ordering 1106 */ 1107 if (a->bs == 4) { 1108 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1109 } 1110 1111 PetscFunctionReturn(0); 1112 } 1113 #undef __FUNC__ 1114 #define __FUNC__ "MatPrintHelp_SeqBAIJ" 1115 int MatPrintHelp_SeqBAIJ(Mat A) 1116 { 1117 static int called = 0; 1118 MPI_Comm comm = A->comm; 1119 1120 PetscFunctionBegin; 1121 if (called) {PetscFunctionReturn(0);} else called = 1; 1122 (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n"); 1123 (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n"); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 /* -------------------------------------------------------------------*/ 1128 static struct _MatOps MatOps = {MatSetValues_SeqBAIJ, 1129 MatGetRow_SeqBAIJ,MatRestoreRow_SeqBAIJ, 1130 MatMult_SeqBAIJ_N,MatMultAdd_SeqBAIJ_N, 1131 MatMultTrans_SeqBAIJ,MatMultTransAdd_SeqBAIJ, 1132 MatSolve_SeqBAIJ_N,0, 1133 0,0, 1134 MatLUFactor_SeqBAIJ,0, 1135 0, 1136 MatTranspose_SeqBAIJ, 1137 MatGetInfo_SeqBAIJ,MatEqual_SeqBAIJ, 1138 MatGetDiagonal_SeqBAIJ,MatDiagonalScale_SeqBAIJ,MatNorm_SeqBAIJ, 1139 0,MatAssemblyEnd_SeqBAIJ, 1140 0, 1141 MatSetOption_SeqBAIJ,MatZeroEntries_SeqBAIJ,MatZeroRows_SeqBAIJ, 1142 MatLUFactorSymbolic_SeqBAIJ,MatLUFactorNumeric_SeqBAIJ_N,0,0, 1143 MatGetSize_SeqBAIJ,MatGetSize_SeqBAIJ,MatGetOwnershipRange_SeqBAIJ, 1144 MatILUFactorSymbolic_SeqBAIJ,0, 1145 0,0, 1146 MatConvertSameType_SeqBAIJ,0,0, 1147 MatILUFactor_SeqBAIJ,0,0, 1148 MatGetSubMatrices_SeqBAIJ,MatIncreaseOverlap_SeqBAIJ, 1149 MatGetValues_SeqBAIJ,0, 1150 MatPrintHelp_SeqBAIJ,MatScale_SeqBAIJ, 1151 0,0,0,MatGetBlockSize_SeqBAIJ, 1152 MatGetRowIJ_SeqBAIJ, 1153 MatRestoreRowIJ_SeqBAIJ, 1154 0,0,0,0,0,0, 1155 MatSetValuesBlocked_SeqBAIJ, 1156 MatGetSubMatrix_SeqBAIJ}; 1157 1158 #undef __FUNC__ 1159 #define __FUNC__ "MatCreateSeqBAIJ" 1160 /*@C 1161 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1162 compressed row) format. For good matrix assembly performance the 1163 user should preallocate the matrix storage by setting the parameter nz 1164 (or the array nzz). By setting these parameters accurately, performance 1165 during matrix assembly can be increased by more than a factor of 50. 1166 1167 Input Parameters: 1168 . comm - MPI communicator, set to PETSC_COMM_SELF 1169 . bs - size of block 1170 . m - number of rows 1171 . n - number of columns 1172 . nz - number of block nonzeros per block row (same for all rows) 1173 . nzz - array containing the number of block nonzeros in the various block rows 1174 (possibly different for each block row) or PETSC_NULL 1175 1176 Output Parameter: 1177 . A - the matrix 1178 1179 Options Database Keys: 1180 $ -mat_no_unroll - uses code that does not unroll the loops in the 1181 $ block calculations (much slower) 1182 $ -mat_block_size - size of the blocks to use 1183 1184 Notes: 1185 The block AIJ format is fully compatible with standard Fortran 77 1186 storage. That is, the stored row and column indices can begin at 1187 either one (as in Fortran) or zero. See the users' manual for details. 1188 1189 Specify the preallocated storage with either nz or nnz (not both). 1190 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1191 allocation. For additional details, see the users manual chapter on 1192 matrices. 1193 1194 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues() 1195 @*/ 1196 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz, Mat *A) 1197 { 1198 Mat B; 1199 Mat_SeqBAIJ *b; 1200 int i,len,ierr,flg,mbs=m/bs,nbs=n/bs,bs2=bs*bs,size; 1201 1202 PetscFunctionBegin; 1203 MPI_Comm_size(comm,&size); 1204 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1205 1206 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,&flg);CHKERRQ(ierr); 1207 1208 if (mbs*bs!=m || nbs*bs!=n) { 1209 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1210 } 1211 1212 *A = 0; 1213 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,comm,MatDestroy,MatView); 1214 PLogObjectCreate(B); 1215 B->data = (void *) (b = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(b); 1216 PetscMemzero(b,sizeof(Mat_SeqBAIJ)); 1217 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1218 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg); CHKERRQ(ierr); 1219 if (!flg) { 1220 switch (bs) { 1221 case 1: 1222 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1223 B->ops->solve = MatSolve_SeqBAIJ_1; 1224 B->ops->mult = MatMult_SeqBAIJ_1; 1225 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1226 break; 1227 case 2: 1228 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1229 B->ops->solve = MatSolve_SeqBAIJ_2; 1230 B->ops->mult = MatMult_SeqBAIJ_2; 1231 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1232 break; 1233 case 3: 1234 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1235 B->ops->solve = MatSolve_SeqBAIJ_3; 1236 B->ops->mult = MatMult_SeqBAIJ_3; 1237 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1238 break; 1239 case 4: 1240 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1241 B->ops->solve = MatSolve_SeqBAIJ_4; 1242 B->ops->mult = MatMult_SeqBAIJ_4; 1243 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1244 break; 1245 case 5: 1246 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1247 B->ops->solve = MatSolve_SeqBAIJ_5; 1248 B->ops->mult = MatMult_SeqBAIJ_5; 1249 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1250 break; 1251 case 7: 1252 B->ops->mult = MatMult_SeqBAIJ_7; 1253 B->ops->solve = MatSolve_SeqBAIJ_7; 1254 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1255 break; 1256 } 1257 } 1258 B->destroy = MatDestroy_SeqBAIJ; 1259 B->view = MatView_SeqBAIJ; 1260 B->factor = 0; 1261 B->lupivotthreshold = 1.0; 1262 B->mapping = 0; 1263 b->row = 0; 1264 b->col = 0; 1265 b->icol = 0; 1266 b->reallocs = 0; 1267 1268 b->m = m; B->m = m; B->M = m; 1269 b->n = n; B->n = n; B->N = n; 1270 b->mbs = mbs; 1271 b->nbs = nbs; 1272 b->imax = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(b->imax); 1273 if (nnz == PETSC_NULL) { 1274 if (nz == PETSC_DEFAULT) nz = 5; 1275 else if (nz <= 0) nz = 1; 1276 for ( i=0; i<mbs; i++ ) b->imax[i] = nz; 1277 nz = nz*mbs; 1278 } else { 1279 nz = 0; 1280 for ( i=0; i<mbs; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1281 } 1282 1283 /* allocate the matrix space */ 1284 len = nz*sizeof(int) + nz*bs2*sizeof(Scalar) + (b->m+1)*sizeof(int); 1285 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1286 PetscMemzero(b->a,nz*bs2*sizeof(Scalar)); 1287 b->j = (int *) (b->a + nz*bs2); 1288 PetscMemzero(b->j,nz*sizeof(int)); 1289 b->i = b->j + nz; 1290 b->singlemalloc = 1; 1291 1292 b->i[0] = 0; 1293 for (i=1; i<mbs+1; i++) { 1294 b->i[i] = b->i[i-1] + b->imax[i-1]; 1295 } 1296 1297 /* b->ilen will count nonzeros in each block row so far. */ 1298 b->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); 1299 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1300 for ( i=0; i<mbs; i++ ) { b->ilen[i] = 0;} 1301 1302 b->bs = bs; 1303 b->bs2 = bs2; 1304 b->mbs = mbs; 1305 b->nz = 0; 1306 b->maxnz = nz*bs2; 1307 b->sorted = 0; 1308 b->roworiented = 1; 1309 b->nonew = 0; 1310 b->diag = 0; 1311 b->solve_work = 0; 1312 b->mult_work = 0; 1313 b->spptr = 0; 1314 B->info.nz_unneeded = (double)b->maxnz; 1315 1316 *A = B; 1317 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1318 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1319 PetscFunctionReturn(0); 1320 } 1321 1322 #undef __FUNC__ 1323 #define __FUNC__ "MatConvertSameType_SeqBAIJ" 1324 int MatConvertSameType_SeqBAIJ(Mat A,Mat *B,int cpvalues) 1325 { 1326 Mat C; 1327 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ *) A->data; 1328 int i,len, mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1329 1330 PetscFunctionBegin; 1331 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1332 1333 *B = 0; 1334 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,A->comm,MatDestroy,MatView); 1335 PLogObjectCreate(C); 1336 C->data = (void *) (c = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(c); 1337 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1338 C->destroy = MatDestroy_SeqBAIJ; 1339 C->view = MatView_SeqBAIJ; 1340 C->factor = A->factor; 1341 c->row = 0; 1342 c->col = 0; 1343 C->assembled = PETSC_TRUE; 1344 1345 c->m = C->m = a->m; 1346 c->n = C->n = a->n; 1347 C->M = a->m; 1348 C->N = a->n; 1349 1350 c->bs = a->bs; 1351 c->bs2 = a->bs2; 1352 c->mbs = a->mbs; 1353 c->nbs = a->nbs; 1354 1355 c->imax = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->imax); 1356 c->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->ilen); 1357 for ( i=0; i<mbs; i++ ) { 1358 c->imax[i] = a->imax[i]; 1359 c->ilen[i] = a->ilen[i]; 1360 } 1361 1362 /* allocate the matrix space */ 1363 c->singlemalloc = 1; 1364 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(Scalar) + sizeof(int)); 1365 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1366 c->j = (int *) (c->a + nz*bs2); 1367 c->i = c->j + nz; 1368 PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int)); 1369 if (mbs > 0) { 1370 PetscMemcpy(c->j,a->j,nz*sizeof(int)); 1371 if (cpvalues == COPY_VALUES) { 1372 PetscMemcpy(c->a,a->a,bs2*nz*sizeof(Scalar)); 1373 } 1374 } 1375 1376 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1377 c->sorted = a->sorted; 1378 c->roworiented = a->roworiented; 1379 c->nonew = a->nonew; 1380 1381 if (a->diag) { 1382 c->diag = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(c->diag); 1383 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1384 for ( i=0; i<mbs; i++ ) { 1385 c->diag[i] = a->diag[i]; 1386 } 1387 } 1388 else c->diag = 0; 1389 c->nz = a->nz; 1390 c->maxnz = a->maxnz; 1391 c->solve_work = 0; 1392 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1393 c->mult_work = 0; 1394 *B = C; 1395 PetscFunctionReturn(0); 1396 } 1397 1398 #undef __FUNC__ 1399 #define __FUNC__ "MatLoad_SeqBAIJ" 1400 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1401 { 1402 Mat_SeqBAIJ *a; 1403 Mat B; 1404 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1,flg; 1405 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1406 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1407 int *masked, nmask,tmp,bs2,ishift; 1408 Scalar *aa; 1409 MPI_Comm comm = ((PetscObject) viewer)->comm; 1410 1411 PetscFunctionBegin; 1412 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,&flg);CHKERRQ(ierr); 1413 bs2 = bs*bs; 1414 1415 MPI_Comm_size(comm,&size); 1416 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1417 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1418 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1419 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1420 M = header[1]; N = header[2]; nz = header[3]; 1421 1422 if (header[3] < 0) { 1423 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1424 } 1425 1426 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1427 1428 /* 1429 This code adds extra rows to make sure the number of rows is 1430 divisible by the blocksize 1431 */ 1432 mbs = M/bs; 1433 extra_rows = bs - M + bs*(mbs); 1434 if (extra_rows == bs) extra_rows = 0; 1435 else mbs++; 1436 if (extra_rows) { 1437 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1438 } 1439 1440 /* read in row lengths */ 1441 rowlengths = (int*) PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1442 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1443 for ( i=0; i<extra_rows; i++ ) rowlengths[M+i] = 1; 1444 1445 /* read in column indices */ 1446 jj = (int*) PetscMalloc( (nz+extra_rows)*sizeof(int) ); CHKPTRQ(jj); 1447 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT); CHKERRQ(ierr); 1448 for ( i=0; i<extra_rows; i++ ) jj[nz+i] = M+i; 1449 1450 /* loop over row lengths determining block row lengths */ 1451 browlengths = (int *) PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1452 PetscMemzero(browlengths,mbs*sizeof(int)); 1453 mask = (int *) PetscMalloc( 2*mbs*sizeof(int) ); CHKPTRQ(mask); 1454 PetscMemzero(mask,mbs*sizeof(int)); 1455 masked = mask + mbs; 1456 rowcount = 0; nzcount = 0; 1457 for ( i=0; i<mbs; i++ ) { 1458 nmask = 0; 1459 for ( j=0; j<bs; j++ ) { 1460 kmax = rowlengths[rowcount]; 1461 for ( k=0; k<kmax; k++ ) { 1462 tmp = jj[nzcount++]/bs; 1463 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1464 } 1465 rowcount++; 1466 } 1467 browlengths[i] += nmask; 1468 /* zero out the mask elements we set */ 1469 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1470 } 1471 1472 /* create our matrix */ 1473 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1474 B = *A; 1475 a = (Mat_SeqBAIJ *) B->data; 1476 1477 /* set matrix "i" values */ 1478 a->i[0] = 0; 1479 for ( i=1; i<= mbs; i++ ) { 1480 a->i[i] = a->i[i-1] + browlengths[i-1]; 1481 a->ilen[i-1] = browlengths[i-1]; 1482 } 1483 a->nz = 0; 1484 for ( i=0; i<mbs; i++ ) a->nz += browlengths[i]; 1485 1486 /* read in nonzero values */ 1487 aa = (Scalar *) PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1488 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR); CHKERRQ(ierr); 1489 for ( i=0; i<extra_rows; i++ ) aa[nz+i] = 1.0; 1490 1491 /* set "a" and "j" values into matrix */ 1492 nzcount = 0; jcount = 0; 1493 for ( i=0; i<mbs; i++ ) { 1494 nzcountb = nzcount; 1495 nmask = 0; 1496 for ( j=0; j<bs; j++ ) { 1497 kmax = rowlengths[i*bs+j]; 1498 for ( k=0; k<kmax; k++ ) { 1499 tmp = jj[nzcount++]/bs; 1500 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1501 } 1502 rowcount++; 1503 } 1504 /* sort the masked values */ 1505 PetscSortInt(nmask,masked); 1506 1507 /* set "j" values into matrix */ 1508 maskcount = 1; 1509 for ( j=0; j<nmask; j++ ) { 1510 a->j[jcount++] = masked[j]; 1511 mask[masked[j]] = maskcount++; 1512 } 1513 /* set "a" values into matrix */ 1514 ishift = bs2*a->i[i]; 1515 for ( j=0; j<bs; j++ ) { 1516 kmax = rowlengths[i*bs+j]; 1517 for ( k=0; k<kmax; k++ ) { 1518 tmp = jj[nzcountb]/bs ; 1519 block = mask[tmp] - 1; 1520 point = jj[nzcountb] - bs*tmp; 1521 idx = ishift + bs2*block + j + bs*point; 1522 a->a[idx] = aa[nzcountb++]; 1523 } 1524 } 1525 /* zero out the mask elements we set */ 1526 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1527 } 1528 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1529 1530 PetscFree(rowlengths); 1531 PetscFree(browlengths); 1532 PetscFree(aa); 1533 PetscFree(jj); 1534 PetscFree(mask); 1535 1536 B->assembled = PETSC_TRUE; 1537 1538 ierr = MatView_Private(B); CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 1543 1544