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