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