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