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