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