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