1 /*$Id: sbaij.c,v 1.1 2000/06/21 15:47:00 balay Exp balay $*/ 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 v_i = *v + M*bs; 273 cols_i = *cols + M*bs; 274 for (i=0; i<bn; i++){ /* for each block row */ 275 M = ai[i+1] - ai[i]; 276 for (j=0; j<M; j++){ 277 itmp = aj[ai[i] + j]; /* block column value */ 278 if (itmp == bn){ 279 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 280 for (k=0; k<bs; k++) { 281 *cols_i++ = i*bs+k; 282 *v_i++ = aa_i[k]; 283 } 284 *ncols += bs; 285 break; 286 } 287 } 288 } 289 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNC__ 294 #define __FUNC__ "MatRestoreRow_SeqSBAIJ" 295 int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 296 { 297 int ierr; 298 299 PetscFunctionBegin; 300 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 301 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 302 PetscFunctionReturn(0); 303 } 304 305 #undef __FUNC__ 306 #define __FUNC__ "MatTranspose_SeqSBAIJ" 307 int MatTranspose_SeqSBAIJ(Mat A,Mat *B) 308 { 309 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 310 Mat C; 311 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 312 int *rows,*cols,bs2=a->bs2,refct; 313 MatScalar *array = a->a; 314 315 PetscFunctionBegin; 316 if (!B && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 317 col = (int*)PetscMalloc((1+nbs)*sizeof(int));CHKPTRQ(col); 318 ierr = PetscMemzero(col,(1+nbs)*sizeof(int));CHKERRQ(ierr); 319 320 for (i=0; i<ai[mbs]; i++) col[aj[i]] += 1; 321 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C);CHKERRQ(ierr); 322 ierr = PetscFree(col);CHKERRQ(ierr); 323 rows = (int*)PetscMalloc(2*bs*sizeof(int));CHKPTRQ(rows); 324 cols = rows + bs; 325 for (i=0; i<mbs; i++) { 326 cols[0] = i*bs; 327 for (k=1; k<bs; k++) cols[k] = cols[k-1] + 1; 328 len = ai[i+1] - ai[i]; 329 for (j=0; j<len; j++) { 330 rows[0] = (*aj++)*bs; 331 for (k=1; k<bs; k++) rows[k] = rows[k-1] + 1; 332 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES);CHKERRQ(ierr); 333 array += bs2; 334 } 335 } 336 ierr = PetscFree(rows);CHKERRQ(ierr); 337 338 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 339 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 340 341 if (B) { 342 *B = C; 343 } else { 344 PetscOps *Abops; 345 MatOps Aops; 346 347 /* This isn't really an in-place transpose */ 348 ierr = PetscFree(a->a);CHKERRQ(ierr); 349 if (!a->singlemalloc) { 350 ierr = PetscFree(a->i);CHKERRQ(ierr); 351 ierr = PetscFree(a->j);CHKERRQ(ierr); 352 } 353 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 354 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 355 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 356 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 357 ierr = PetscFree(a);CHKERRQ(ierr); 358 359 360 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 361 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 362 363 /* 364 This is horrible, horrible code. We need to keep the 365 A pointers for the bops and ops but copy everything 366 else from C. 367 */ 368 Abops = A->bops; 369 Aops = A->ops; 370 refct = A->refct; 371 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 372 A->bops = Abops; 373 A->ops = Aops; 374 A->qlist = 0; 375 A->refct = refct; 376 377 PetscHeaderDestroy(C); 378 } 379 PetscFunctionReturn(0); 380 } 381 382 #undef __FUNC__ 383 #define __FUNC__ "MatView_SeqSBAIJ_Binary" 384 static int MatView_SeqSBAIJ_Binary(Mat A,Viewer viewer) 385 { 386 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 387 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 388 Scalar *aa; 389 FILE *file; 390 391 PetscFunctionBegin; 392 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 393 col_lens = (int*)PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 394 col_lens[0] = MAT_COOKIE; 395 396 col_lens[1] = a->m; 397 col_lens[2] = a->m; 398 col_lens[3] = a->s_nz*bs2; 399 400 /* store lengths of each row and write (including header) to file */ 401 count = 0; 402 for (i=0; i<a->mbs; i++) { 403 for (j=0; j<bs; j++) { 404 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 405 } 406 } 407 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1);CHKERRQ(ierr); 408 ierr = PetscFree(col_lens);CHKERRQ(ierr); 409 410 /* store column indices (zero start index) */ 411 jj = (int*)PetscMalloc((a->s_nz+1)*bs2*sizeof(int));CHKPTRQ(jj); 412 count = 0; 413 for (i=0; i<a->mbs; i++) { 414 for (j=0; j<bs; j++) { 415 for (k=a->i[i]; k<a->i[i+1]; k++) { 416 for (l=0; l<bs; l++) { 417 jj[count++] = bs*a->j[k] + l; 418 } 419 } 420 } 421 } 422 ierr = PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);CHKERRQ(ierr); 423 ierr = PetscFree(jj);CHKERRQ(ierr); 424 425 /* store nonzero values */ 426 aa = (Scalar*)PetscMalloc((a->s_nz+1)*bs2*sizeof(Scalar));CHKPTRQ(aa); 427 count = 0; 428 for (i=0; i<a->mbs; i++) { 429 for (j=0; j<bs; j++) { 430 for (k=a->i[i]; k<a->i[i+1]; k++) { 431 for (l=0; l<bs; l++) { 432 aa[count++] = a->a[bs2*k + l*bs + j]; 433 } 434 } 435 } 436 } 437 ierr = PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);CHKERRQ(ierr); 438 ierr = PetscFree(aa);CHKERRQ(ierr); 439 440 ierr = ViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 441 if (file) { 442 fprintf(file,"-matload_block_size %d\n",a->bs); 443 } 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNC__ 448 #define __FUNC__ "MatView_SeqSBAIJ_ASCII" 449 static int MatView_SeqSBAIJ_ASCII(Mat A,Viewer viewer) 450 { 451 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 452 int ierr,i,j,format,bs = a->bs,k,l,bs2=a->bs2; 453 char *outputname; 454 455 PetscFunctionBegin; 456 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 457 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 458 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 459 ierr = ViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 460 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 461 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 462 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 463 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 464 for (i=0; i<a->mbs; i++) { 465 for (j=0; j<bs; j++) { 466 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 467 for (k=a->i[i]; k<a->i[i+1]; k++) { 468 for (l=0; l<bs; l++) { 469 #if defined(PETSC_USE_COMPLEX) 470 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 471 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 472 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 473 } else 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 (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 477 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 478 } 479 #else 480 if (a->a[bs2*k + l*bs + j] != 0.0) { 481 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 482 } 483 #endif 484 } 485 } 486 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 487 } 488 } 489 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 490 } else { 491 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 492 for (i=0; i<a->mbs; i++) { 493 for (j=0; j<bs; j++) { 494 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 495 for (k=a->i[i]; k<a->i[i+1]; k++) { 496 for (l=0; l<bs; l++) { 497 #if defined(PETSC_USE_COMPLEX) 498 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 499 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 500 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 501 } else 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 { 505 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 506 } 507 #else 508 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 509 #endif 510 } 511 } 512 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 513 } 514 } 515 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 516 } 517 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNC__ 522 #define __FUNC__ "MatView_SeqSBAIJ_Draw_Zoom" 523 static int MatView_SeqSBAIJ_Draw_Zoom(Draw draw,void *Aa) 524 { 525 Mat A = (Mat) Aa; 526 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 527 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 528 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 529 MatScalar *aa; 530 MPI_Comm comm; 531 Viewer viewer; 532 533 PetscFunctionBegin; 534 /* 535 This is nasty. If this is called from an originally parallel matrix 536 then all processes call this,but only the first has the matrix so the 537 rest should return immediately. 538 */ 539 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 540 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 541 if (rank) PetscFunctionReturn(0); 542 543 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 544 545 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 546 DrawString(draw, .3*(xl+xr), .3*(yl+yr), DRAW_BLACK, "symmetric"); 547 548 /* loop over matrix elements drawing boxes */ 549 color = DRAW_BLUE; 550 for (i=0,row=0; i<mbs; i++,row+=bs) { 551 for (j=a->i[i]; j<a->i[i+1]; j++) { 552 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 553 x_l = a->j[j]*bs; x_r = x_l + 1.0; 554 aa = a->a + j*bs2; 555 for (k=0; k<bs; k++) { 556 for (l=0; l<bs; l++) { 557 if (PetscRealPart(*aa++) >= 0.) continue; 558 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 559 } 560 } 561 } 562 } 563 color = DRAW_CYAN; 564 for (i=0,row=0; i<mbs; i++,row+=bs) { 565 for (j=a->i[i]; j<a->i[i+1]; j++) { 566 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 567 x_l = a->j[j]*bs; x_r = x_l + 1.0; 568 aa = a->a + j*bs2; 569 for (k=0; k<bs; k++) { 570 for (l=0; l<bs; l++) { 571 if (PetscRealPart(*aa++) != 0.) continue; 572 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 573 } 574 } 575 } 576 } 577 578 color = DRAW_RED; 579 for (i=0,row=0; i<mbs; i++,row+=bs) { 580 for (j=a->i[i]; j<a->i[i+1]; j++) { 581 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 582 x_l = a->j[j]*bs; x_r = x_l + 1.0; 583 aa = a->a + j*bs2; 584 for (k=0; k<bs; k++) { 585 for (l=0; l<bs; l++) { 586 if (PetscRealPart(*aa++) <= 0.) continue; 587 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 588 } 589 } 590 } 591 } 592 PetscFunctionReturn(0); 593 } 594 595 #undef __FUNC__ 596 #define __FUNC__ "MatView_SeqSBAIJ_Draw" 597 static int MatView_SeqSBAIJ_Draw(Mat A,Viewer viewer) 598 { 599 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 600 int ierr; 601 PetscReal xl,yl,xr,yr,w,h; 602 Draw draw; 603 PetscTruth isnull; 604 605 PetscFunctionBegin; 606 607 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 608 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 609 610 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 611 xr = a->m; yr = a->m; h = yr/10.0; w = xr/10.0; 612 xr += w; yr += h; xl = -w; yl = -h; 613 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 614 ierr = DrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 615 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 616 PetscFunctionReturn(0); 617 } 618 619 #undef __FUNC__ 620 #define __FUNC__ "MatView_SeqSBAIJ" 621 int MatView_SeqSBAIJ(Mat A,Viewer viewer) 622 { 623 int ierr; 624 PetscTruth issocket,isascii,isbinary,isdraw; 625 626 PetscFunctionBegin; 627 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 628 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 629 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 630 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 631 if (issocket) { 632 SETERRQ(PETSC_ERR_SUP,0,"Socket viewer not supported"); 633 } else if (isascii){ 634 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 635 } else if (isbinary) { 636 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 637 } else if (isdraw) { 638 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 639 } else { 640 SETERRQ1(1,1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 641 } 642 PetscFunctionReturn(0); 643 } 644 645 646 #undef __FUNC__ 647 #define __FUNC__ "MatGetValues_SeqSBAIJ" 648 int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 649 { 650 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 651 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 652 int *ai = a->i,*ailen = a->ilen; 653 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 654 MatScalar *ap,*aa = a->a,zero = 0.0; 655 656 PetscFunctionBegin; 657 for (k=0; k<m; k++) { /* loop over rows */ 658 row = im[k]; brow = row/bs; 659 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 660 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 661 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 662 nrow = ailen[brow]; 663 for (l=0; l<n; l++) { /* loop over columns */ 664 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 665 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 666 col = in[l] ; 667 bcol = col/bs; 668 cidx = col%bs; 669 ridx = row%bs; 670 high = nrow; 671 low = 0; /* assume unsorted */ 672 while (high-low > 5) { 673 t = (low+high)/2; 674 if (rp[t] > bcol) high = t; 675 else low = t; 676 } 677 for (i=low; i<high; i++) { 678 if (rp[i] > bcol) break; 679 if (rp[i] == bcol) { 680 *v++ = ap[bs2*i+bs*cidx+ridx]; 681 goto finished; 682 } 683 } 684 *v++ = zero; 685 finished:; 686 } 687 } 688 PetscFunctionReturn(0); 689 } 690 691 692 #undef __FUNC__ 693 #define __FUNC__ "MatSetValuesBlocked_SeqSBAIJ" 694 int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 695 { 696 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 697 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 698 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 699 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval,ierr; 700 Scalar *value = v; 701 MatScalar *ap,*aa=a->a,*bap; 702 703 PetscFunctionBegin; 704 if (roworiented) { 705 stepval = (n-1)*bs; 706 } else { 707 stepval = (m-1)*bs; 708 } 709 for (k=0; k<m; k++) { /* loop over added rows */ 710 row = im[k]; 711 if (row < 0) continue; 712 #if defined(PETSC_USE_BOPT_g) 713 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 714 #endif 715 rp = aj + ai[row]; 716 ap = aa + bs2*ai[row]; 717 rmax = imax[row]; 718 nrow = ailen[row]; 719 low = 0; 720 for (l=0; l<n; l++) { /* loop over added columns */ 721 if (in[l] < 0) continue; 722 #if defined(PETSC_USE_BOPT_g) 723 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 724 #endif 725 col = in[l]; 726 if (roworiented) { 727 value = v + k*(stepval+bs)*bs + l*bs; 728 } else { 729 value = v + l*(stepval+bs)*bs + k*bs; 730 } 731 if (!sorted) low = 0; high = nrow; 732 while (high-low > 7) { 733 t = (low+high)/2; 734 if (rp[t] > col) high = t; 735 else low = t; 736 } 737 for (i=low; i<high; i++) { 738 if (rp[i] > col) break; 739 if (rp[i] == col) { 740 bap = ap + bs2*i; 741 if (roworiented) { 742 if (is == ADD_VALUES) { 743 for (ii=0; ii<bs; ii++,value+=stepval) { 744 for (jj=ii; jj<bs2; jj+=bs) { 745 bap[jj] += *value++; 746 } 747 } 748 } else { 749 for (ii=0; ii<bs; ii++,value+=stepval) { 750 for (jj=ii; jj<bs2; jj+=bs) { 751 bap[jj] = *value++; 752 } 753 } 754 } 755 } else { 756 if (is == ADD_VALUES) { 757 for (ii=0; ii<bs; ii++,value+=stepval) { 758 for (jj=0; jj<bs; jj++) { 759 *bap++ += *value++; 760 } 761 } 762 } else { 763 for (ii=0; ii<bs; ii++,value+=stepval) { 764 for (jj=0; jj<bs; jj++) { 765 *bap++ = *value++; 766 } 767 } 768 } 769 } 770 goto noinsert2; 771 } 772 } 773 if (nonew == 1) goto noinsert2; 774 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 775 if (nrow >= rmax) { 776 /* there is no extra room in row, therefore enlarge */ 777 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 778 MatScalar *new_a; 779 780 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 781 782 /* malloc new storage space */ 783 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 784 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 785 new_j = (int*)(new_a + bs2*new_nz); 786 new_i = new_j + new_nz; 787 788 /* copy over old data into new slots */ 789 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 790 for (ii=row+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 791 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 792 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 793 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 794 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 795 ierr = PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 796 ierr = PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 797 /* free up old matrix storage */ 798 ierr = PetscFree(a->a);CHKERRQ(ierr); 799 if (!a->singlemalloc) { 800 ierr = PetscFree(a->i);CHKERRQ(ierr); 801 ierr = PetscFree(a->j);CHKERRQ(ierr); 802 } 803 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 804 a->singlemalloc = PETSC_TRUE; 805 806 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 807 rmax = imax[row] = imax[row] + CHUNKSIZE; 808 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 809 a->maxnz += bs2*CHUNKSIZE; 810 a->reallocs++; 811 a->nz++; 812 } 813 N = nrow++ - 1; 814 /* shift up all the later entries in this row */ 815 for (ii=N; ii>=i; ii--) { 816 rp[ii+1] = rp[ii]; 817 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 818 } 819 if (N >= i) { 820 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 821 } 822 rp[i] = col; 823 bap = ap + bs2*i; 824 if (roworiented) { 825 for (ii=0; ii<bs; ii++,value+=stepval) { 826 for (jj=ii; jj<bs2; jj+=bs) { 827 bap[jj] = *value++; 828 } 829 } 830 } else { 831 for (ii=0; ii<bs; ii++,value+=stepval) { 832 for (jj=0; jj<bs; jj++) { 833 *bap++ = *value++; 834 } 835 } 836 } 837 noinsert2:; 838 low = i; 839 } 840 ailen[row] = nrow; 841 } 842 PetscFunctionReturn(0); 843 } 844 845 #undef __FUNC__ 846 #define __FUNC__ "MatAssemblyEnd_SeqSBAIJ" 847 int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 848 { 849 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 850 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 851 int m = a->m,*ip,N,*ailen = a->ilen; 852 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 853 MatScalar *aa = a->a,*ap; 854 855 PetscFunctionBegin; 856 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 857 858 if (m) rmax = ailen[0]; 859 for (i=1; i<mbs; i++) { 860 /* move each row back by the amount of empty slots (fshift) before it*/ 861 fshift += imax[i-1] - ailen[i-1]; 862 rmax = PetscMax(rmax,ailen[i]); 863 if (fshift) { 864 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 865 N = ailen[i]; 866 for (j=0; j<N; j++) { 867 ip[j-fshift] = ip[j]; 868 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 869 } 870 } 871 ai[i] = ai[i-1] + ailen[i-1]; 872 } 873 if (mbs) { 874 fshift += imax[mbs-1] - ailen[mbs-1]; 875 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 876 } 877 /* reset ilen and imax for each row */ 878 for (i=0; i<mbs; i++) { 879 ailen[i] = imax[i] = ai[i+1] - ai[i]; 880 } 881 a->s_nz = ai[mbs]; 882 883 /* diagonals may have moved, so kill the diagonal pointers */ 884 if (fshift && a->diag) { 885 ierr = PetscFree(a->diag);CHKERRQ(ierr); 886 PLogObjectMemory(A,-(m+1)*sizeof(int)); 887 a->diag = 0; 888 } 889 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 890 m,a->m,a->bs,fshift*bs2,a->s_nz*bs2); 891 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 892 a->reallocs); 893 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 894 a->reallocs = 0; 895 A->info.nz_unneeded = (PetscReal)fshift*bs2; 896 897 PetscFunctionReturn(0); 898 } 899 900 /* 901 This function returns an array of flags which indicate the locations of contiguous 902 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 903 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 904 Assume: sizes should be long enough to hold all the values. 905 */ 906 #undef __FUNC__ 907 #define __FUNC__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 908 static int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 909 { 910 int i,j,k,row; 911 PetscTruth flg; 912 913 PetscFunctionBegin; 914 for (i=0,j=0; i<n; j++) { 915 row = idx[i]; 916 if (row%bs!=0) { /* Not the begining of a block */ 917 sizes[j] = 1; 918 i++; 919 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 920 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 921 i++; 922 } else { /* Begining of the block, so check if the complete block exists */ 923 flg = PETSC_TRUE; 924 for (k=1; k<bs; k++) { 925 if (row+k != idx[i+k]) { /* break in the block */ 926 flg = PETSC_FALSE; 927 break; 928 } 929 } 930 if (flg == PETSC_TRUE) { /* No break in the bs */ 931 sizes[j] = bs; 932 i+= bs; 933 } else { 934 sizes[j] = 1; 935 i++; 936 } 937 } 938 } 939 *bs_max = j; 940 PetscFunctionReturn(0); 941 } 942 943 #undef __FUNC__ 944 #define __FUNC__ "MatZeroRows_SeqSBAIJ" 945 int MatZeroRows_SeqSBAIJ(Mat A,IS is,Scalar *diag) 946 { 947 Mat_SeqSBAIJ *sbaij=(Mat_SeqSBAIJ*)A->data; 948 int ierr,i,j,k,count,is_n,*is_idx,*rows; 949 int bs=sbaij->bs,bs2=sbaij->bs2,*sizes,row,bs_max; 950 Scalar zero = 0.0; 951 MatScalar *aa; 952 953 PetscFunctionBegin; 954 /* Make a copy of the IS and sort it */ 955 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 956 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 957 958 /* allocate memory for rows,sizes */ 959 rows = (int*)PetscMalloc((3*is_n+1)*sizeof(int));CHKPTRQ(rows); 960 sizes = rows + is_n; 961 962 /* initialize copy IS values to rows, and sort them */ 963 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 964 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 965 if (sbaij->keepzeroedrows) { /* do not change nonzero structure */ 966 for (i=0; i<is_n; i++) { sizes[i] = 1; } /* sizes: size of blocks, = 1 or bs */ 967 bs_max = is_n; /* bs_max: num. of contiguous block row in the row */ 968 } else { 969 ierr = MatZeroRows_SeqSBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 970 } 971 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 972 973 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 974 row = rows[j]; /* row to be zeroed */ 975 if (row < 0 || row > sbaij->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"row %d out of range",row); 976 count = (sbaij->i[row/bs +1] - sbaij->i[row/bs])*bs; /* num. of elements in the row */ 977 aa = sbaij->a + sbaij->i[row/bs]*bs2 + (row%bs); 978 if (sizes[i] == bs && !sbaij->keepzeroedrows) { 979 if (diag) { 980 if (sbaij->ilen[row/bs] > 0) { 981 sbaij->ilen[row/bs] = 1; 982 sbaij->j[sbaij->i[row/bs]] = row/bs; 983 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 984 } 985 /* Now insert all the diagoanl values for this bs */ 986 for (k=0; k<bs; k++) { 987 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 988 } 989 } else { /* (!diag) */ 990 sbaij->ilen[row/bs] = 0; 991 } /* end (!diag) */ 992 } else { /* (sizes[i] != bs), broken block */ 993 #if defined (PETSC_USE_DEBUG) 994 if (sizes[i] != 1) SETERRQ(1,0,"Internal Error. Value should be 1"); 995 #endif 996 for (k=0; k<count; k++) { 997 aa[0] = zero; 998 aa+=bs; 999 } 1000 if (diag) { 1001 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 1002 } 1003 } 1004 } 1005 1006 ierr = PetscFree(rows);CHKERRQ(ierr); 1007 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1008 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011 1012 /* Only add/insert a(i,j) with i<=j (blocks). 1013 Any a(i,j) with i>j input by user is ingored. 1014 */ 1015 1016 #undef __FUNC__ 1017 #define __FUNC__ "MatSetValues_SeqSBAIJ" 1018 int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 1019 { 1020 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1021 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 1022 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 1023 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 1024 int ridx,cidx,bs2=a->bs2,ierr; 1025 MatScalar *ap,value,*aa=a->a,*bap; 1026 int ambs=a->mbs; 1027 1028 PetscFunctionBegin; 1029 1030 for (k=0; k<m; k++) { /* loop over added rows */ 1031 row = im[k]; /* row number */ 1032 brow = row/bs; /* block row number */ 1033 if (row < 0) continue; 1034 #if defined(PETSC_USE_BOPT_g) 1035 if (row >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large: row %d max %d",row,a->m); 1036 #endif 1037 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 1038 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 1039 rmax = imax[brow]; /* maximum space allocated for this row */ 1040 nrow = ailen[brow]; /* actual length of this row */ 1041 low = 0; 1042 1043 for (l=0; l<n; l++) { /* loop over added columns */ 1044 if (in[l] < 0) continue; 1045 #if defined(PETSC_USE_BOPT_g) 1046 if (in[l] >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large: col %d max %d",in[l],a->m); 1047 #endif 1048 col = in[l]; 1049 bcol = col/bs; /* block col number */ 1050 1051 if (brow <= bcol){ 1052 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 1053 /* if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ */ 1054 /* element value a(k,l) */ 1055 if (roworiented) { 1056 value = v[l + k*n]; 1057 } else { 1058 value = v[k + l*m]; 1059 } 1060 1061 /* move pointer bap to a(k,l) quickly and add/insert value */ 1062 if (!sorted) low = 0; high = nrow; 1063 while (high-low > 7) { 1064 t = (low+high)/2; 1065 if (rp[t] > bcol) high = t; 1066 else low = t; 1067 } 1068 for (i=low; i<high; i++) { 1069 /* printf("The loop of i=low.., rp[%d]=%d\n",i,rp[i]); */ 1070 if (rp[i] > bcol) break; 1071 if (rp[i] == bcol) { 1072 bap = ap + bs2*i + bs*cidx + ridx; 1073 if (is == ADD_VALUES) *bap += value; 1074 else *bap = value; 1075 goto noinsert1; 1076 } 1077 } 1078 1079 if (nonew == 1) goto noinsert1; 1080 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1081 if (nrow >= rmax) { 1082 /* there is no extra room in row, therefore enlarge */ 1083 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 1084 MatScalar *new_a; 1085 1086 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1087 1088 /* Malloc new storage space */ 1089 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 1090 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 1091 new_j = (int*)(new_a + bs2*new_nz); 1092 new_i = new_j + new_nz; 1093 1094 /* copy over old data into new slots */ 1095 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 1096 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 1097 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 1098 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 1099 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 1100 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1101 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 1102 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 1103 /* free up old matrix storage */ 1104 ierr = PetscFree(a->a);CHKERRQ(ierr); 1105 if (!a->singlemalloc) { 1106 ierr = PetscFree(a->i);CHKERRQ(ierr); 1107 ierr = PetscFree(a->j);CHKERRQ(ierr); 1108 } 1109 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 1110 a->singlemalloc = PETSC_TRUE; 1111 1112 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1113 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1114 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 1115 a->s_maxnz += bs2*CHUNKSIZE; 1116 a->reallocs++; 1117 a->s_nz++; 1118 } 1119 1120 N = nrow++ - 1; 1121 /* shift up all the later entries in this row */ 1122 for (ii=N; ii>=i; ii--) { 1123 rp[ii+1] = rp[ii]; 1124 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1125 } 1126 if (N>=i) { 1127 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1128 } 1129 rp[i] = bcol; 1130 ap[bs2*i + bs*cidx + ridx] = value; 1131 noinsert1:; 1132 low = i; 1133 /* } */ 1134 } /* end of if .. if.. */ 1135 } /* end of loop over added columns */ 1136 ailen[brow] = nrow; 1137 } /* end of loop over added rows */ 1138 1139 PetscFunctionReturn(0); 1140 } 1141 1142 extern int MatLUFactorSymbolic_SeqSBAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1143 extern int MatLUFactor_SeqSBAIJ(Mat,IS,IS,MatLUInfo*); 1144 extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int); 1145 extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 1146 extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 1147 extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec); 1148 extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec); 1149 extern int MatScale_SeqSBAIJ(Scalar*,Mat); 1150 extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 1151 extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 1152 extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec); 1153 extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 1154 extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 1155 extern int MatZeroEntries_SeqSBAIJ(Mat); 1156 1157 extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 1158 extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 1159 extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 1160 extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 1161 extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 1162 extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 1163 extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 1164 extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 1165 extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec); 1166 extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec); 1167 extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec); 1168 extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec); 1169 extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec); 1170 extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec); 1171 extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec); 1172 1173 extern int MatLUFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 1174 extern int MatLUFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 1175 extern int MatLUFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 1176 extern int MatLUFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 1177 extern int MatLUFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 1178 extern int MatLUFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 1179 extern int MatLUFactorNumeric_SeqSBAIJ_6(Mat,Mat*); 1180 1181 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 1182 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 1183 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 1184 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 1185 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 1186 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 1187 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 1188 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 1189 1190 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 1191 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 1192 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 1193 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 1194 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 1195 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 1196 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 1197 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 1198 1199 #undef __FUNC__ 1200 #define __FUNC__ "MatILUFactor_SeqSBAIJ" 1201 int MatILUFactor_SeqSBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1202 { 1203 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)inA->data; 1204 Mat outA; 1205 int ierr; 1206 PetscTruth row_identity,col_identity; 1207 1208 PetscFunctionBegin; 1209 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,0,"Only levels = 0 supported for in-place ILU"); 1210 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1211 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1212 if (!row_identity || !col_identity) { 1213 SETERRQ(1,1,"Row and column permutations must be identity for in-place ILU"); 1214 } 1215 1216 outA = inA; 1217 inA->factor = FACTOR_LU; 1218 1219 if (!a->diag) { 1220 ierr = MatMarkDiagonal_SeqBAIJ(inA);CHKERRQ(ierr); 1221 } 1222 /* 1223 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 1224 for ILU(0) factorization with natural ordering 1225 */ 1226 switch (a->bs) { 1227 case 1: 1228 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1229 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1230 case 2: 1231 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 1232 inA->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 1233 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1234 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1235 break; 1236 case 3: 1237 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 1238 inA->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; 1239 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; 1240 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1241 break; 1242 case 4: 1243 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1244 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1245 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1246 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1247 break; 1248 case 5: 1249 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1250 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1251 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1252 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1253 break; 1254 case 6: 1255 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1256 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1257 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1258 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1259 break; 1260 case 7: 1261 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1262 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1263 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1264 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1265 break; 1266 default: 1267 a->row = row; 1268 a->col = col; 1269 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1270 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1271 1272 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1273 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1274 PLogObjectParent(inA,a->icol); 1275 1276 if (!a->solve_work) { 1277 a->solve_work = (Scalar*)PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1278 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1279 } 1280 } 1281 1282 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1283 1284 PetscFunctionReturn(0); 1285 } 1286 #undef __FUNC__ 1287 #define __FUNC__ "MatPrintHelp_SeqSBAIJ" 1288 int MatPrintHelp_SeqSBAIJ(Mat A) 1289 { 1290 static PetscTruth called = PETSC_FALSE; 1291 MPI_Comm comm = A->comm; 1292 int ierr; 1293 1294 PetscFunctionBegin; 1295 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1296 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1297 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 EXTERN_C_BEGIN 1302 #undef __FUNC__ 1303 #define __FUNC__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1304 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1305 { 1306 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1307 int i,nz,n; 1308 1309 PetscFunctionBegin; 1310 nz = baij->maxnz; 1311 n = baij->n; 1312 for (i=0; i<nz; i++) { 1313 baij->j[i] = indices[i]; 1314 } 1315 baij->nz = nz; 1316 for (i=0; i<n; i++) { 1317 baij->ilen[i] = baij->imax[i]; 1318 } 1319 1320 PetscFunctionReturn(0); 1321 } 1322 EXTERN_C_END 1323 1324 #undef __FUNC__ 1325 #define __FUNC__ "MatSeqSBAIJSetColumnIndices" 1326 /*@ 1327 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1328 in the matrix. 1329 1330 Input Parameters: 1331 + mat - the SeqBAIJ matrix 1332 - indices - the column indices 1333 1334 Level: advanced 1335 1336 Notes: 1337 This can be called if you have precomputed the nonzero structure of the 1338 matrix and want to provide it to the matrix object to improve the performance 1339 of the MatSetValues() operation. 1340 1341 You MUST have set the correct numbers of nonzeros per row in the call to 1342 MatCreateSeqBAIJ(). 1343 1344 MUST be called before any calls to MatSetValues(); 1345 1346 @*/ 1347 int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices) 1348 { 1349 int ierr,(*f)(Mat,int *); 1350 1351 PetscFunctionBegin; 1352 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1353 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1354 if (f) { 1355 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1356 } else { 1357 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1358 } 1359 PetscFunctionReturn(0); 1360 } 1361 1362 /* -------------------------------------------------------------------*/ 1363 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1364 MatGetRow_SeqSBAIJ, 1365 MatRestoreRow_SeqSBAIJ, 1366 MatMult_SeqSBAIJ_N, 1367 MatMultAdd_SeqSBAIJ_N, 1368 MatMultTranspose_SeqSBAIJ, 1369 MatMultTransposeAdd_SeqSBAIJ, 1370 MatSolve_SeqSBAIJ_N, 1371 0, 1372 0, 1373 0, 1374 MatLUFactor_SeqSBAIJ, 1375 0, 1376 0, 1377 MatTranspose_SeqSBAIJ, 1378 MatGetInfo_SeqSBAIJ, 1379 MatEqual_SeqSBAIJ, 1380 MatGetDiagonal_SeqSBAIJ, 1381 MatDiagonalScale_SeqSBAIJ, 1382 MatNorm_SeqSBAIJ, 1383 0, 1384 MatAssemblyEnd_SeqSBAIJ, 1385 0, 1386 MatSetOption_SeqSBAIJ, 1387 MatZeroEntries_SeqSBAIJ, 1388 MatZeroRows_SeqSBAIJ, 1389 MatLUFactorSymbolic_SeqSBAIJ, 1390 MatLUFactorNumeric_SeqSBAIJ_N, 1391 0, 1392 0, 1393 MatGetSize_SeqSBAIJ, 1394 MatGetSize_SeqSBAIJ, 1395 MatGetOwnershipRange_SeqSBAIJ, 1396 MatILUFactorSymbolic_SeqSBAIJ, 1397 0, 1398 0, 1399 0, 1400 MatDuplicate_SeqSBAIJ, 1401 0, 1402 0, 1403 MatILUFactor_SeqSBAIJ, 1404 0, 1405 0, 1406 MatGetSubMatrices_SeqSBAIJ, 1407 MatIncreaseOverlap_SeqSBAIJ, 1408 MatGetValues_SeqSBAIJ, 1409 0, 1410 MatPrintHelp_SeqSBAIJ, 1411 MatScale_SeqSBAIJ, 1412 0, 1413 0, 1414 0, 1415 MatGetBlockSize_SeqSBAIJ, 1416 MatGetRowIJ_SeqSBAIJ, 1417 MatRestoreRowIJ_SeqSBAIJ, 1418 0, 1419 0, 1420 0, 1421 0, 1422 0, 1423 0, 1424 MatSetValuesBlocked_SeqSBAIJ, 1425 MatGetSubMatrix_SeqSBAIJ, 1426 0, 1427 0, 1428 MatGetMaps_Petsc}; 1429 1430 EXTERN_C_BEGIN 1431 #undef __FUNC__ 1432 #define __FUNC__ "MatStoreValues_SeqSBAIJ" 1433 int MatStoreValues_SeqSBAIJ(Mat mat) 1434 { 1435 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1436 int nz = aij->i[aij->m]*aij->bs*aij->bs2; 1437 int ierr; 1438 1439 PetscFunctionBegin; 1440 if (aij->nonew != 1) { 1441 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1442 } 1443 1444 /* allocate space for values if not already there */ 1445 if (!aij->saved_values) { 1446 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1447 } 1448 1449 /* copy values over */ 1450 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 EXTERN_C_END 1454 1455 EXTERN_C_BEGIN 1456 #undef __FUNC__ 1457 #define __FUNC__ "MatRetrieveValues_SeqSBAIJ" 1458 int MatRetrieveValues_SeqSBAIJ(Mat mat) 1459 { 1460 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1461 int nz = aij->i[aij->m]*aij->bs*aij->bs2,ierr; 1462 1463 PetscFunctionBegin; 1464 if (aij->nonew != 1) { 1465 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1466 } 1467 if (!aij->saved_values) { 1468 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1469 } 1470 1471 /* copy values over */ 1472 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1473 PetscFunctionReturn(0); 1474 } 1475 EXTERN_C_END 1476 1477 /*------------------------------------ New 1 -------------------------------*/ 1478 #undef __FUNC__ 1479 #define __FUNC__ "MatCreateSeqSBAIJ" 1480 /*@C 1481 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1482 compressed row) format. For good matrix assembly performance the 1483 user should preallocate the matrix storage by setting the parameter nz 1484 (or the array nnz). By setting these parameters accurately, performance 1485 during matrix assembly can be increased by more than a factor of 50. 1486 1487 Collective on MPI_Comm 1488 1489 Input Parameters: 1490 + comm - MPI communicator, set to PETSC_COMM_SELF 1491 . bs - size of block 1492 . m - number of rows, or number of columns 1493 . nz - number of block nonzeros per block row (same for all rows) 1494 - nnz - array containing the number of block nonzeros in the various block rows 1495 (possibly different for each block row) or PETSC_NULL 1496 1497 Output Parameter: 1498 . A - the symmetric matrix 1499 1500 Options Database Keys: 1501 . -mat_no_unroll - uses code that does not unroll the loops in the 1502 block calculations (much slower) 1503 . -mat_block_size - size of the blocks to use 1504 1505 Level: intermediate 1506 1507 Notes: 1508 The block AIJ format is fully compatible with standard Fortran 77 1509 storage. That is, the stored row and column indices can begin at 1510 either one (as in Fortran) or zero. See the users' manual for details. 1511 1512 Specify the preallocated storage with either nz or nnz (not both). 1513 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1514 allocation. For additional details, see the users manual chapter on 1515 matrices. 1516 1517 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1518 @*/ 1519 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1520 { 1521 Mat B; 1522 Mat_SeqSBAIJ *b; 1523 int i,len,ierr,mbs,nbs,bs2,size; 1524 PetscTruth flg; 1525 int s_nz,s_maxnz,*imax; 1526 1527 PetscFunctionBegin; 1528 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1529 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1530 1531 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1532 mbs = m/bs; 1533 /* nbs = n/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 = nbs; */ 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