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