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