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