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