1 /*$Id: sbaij.c,v 1.9 2000/07/26 14:07:20 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 MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*); 1145 extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal); 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 MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 1176 extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 1177 extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 1178 extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 1179 extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 1180 extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 1181 extern int MatCholeskyFactorNumeric_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 #ifdef MatIncompleteCholeskyFactor 1202 /* This function is modified from MatILUFactor_SeqSBAIJ. 1203 Needs further work! Don't forget to add the function to the matrix table. */ 1204 #undef __FUNC__ 1205 #define __FUNC__ "MatIncompleteCholeskyFactor_SeqSBAIJ" 1206 int MatIncompleteCholeskyFactor_SeqSBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1207 { 1208 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1209 Mat outA; 1210 int ierr; 1211 PetscTruth row_identity,col_identity; 1212 1213 PetscFunctionBegin; 1214 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,0,"Only levels = 0 supported for in-place ILU"); 1215 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1216 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1217 if (!row_identity || !col_identity) { 1218 SETERRQ(1,1,"Row and column permutations must be identity for in-place ILU"); 1219 } 1220 1221 outA = inA; 1222 inA->factor = FACTOR_LU; 1223 1224 if (!a->diag) { 1225 ierr = MatMarkDiagonal_SeqBAIJ(inA);CHKERRQ(ierr); 1226 } 1227 /* 1228 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 1229 for ILU(0) factorization with natural ordering 1230 */ 1231 switch (a->bs) { 1232 case 1: 1233 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1234 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1235 case 2: 1236 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 1237 inA->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 1238 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1239 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1240 break; 1241 case 3: 1242 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 1243 inA->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; 1244 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; 1245 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1246 break; 1247 case 4: 1248 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1249 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1250 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1251 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1252 break; 1253 case 5: 1254 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1255 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1256 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1257 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1258 break; 1259 case 6: 1260 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1261 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1262 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1263 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1264 break; 1265 case 7: 1266 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1267 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1268 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1269 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1270 break; 1271 default: 1272 a->row = row; 1273 a->col = col; 1274 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1275 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1276 1277 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1278 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1279 PLogObjectParent(inA,a->icol); 1280 1281 if (!a->solve_work) { 1282 a->solve_work = (Scalar*)PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1283 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1284 } 1285 } 1286 1287 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1288 1289 PetscFunctionReturn(0); 1290 } 1291 #endif 1292 1293 #undef __FUNC__ 1294 #define __FUNC__ "MatPrintHelp_SeqSBAIJ" 1295 int MatPrintHelp_SeqSBAIJ(Mat A) 1296 { 1297 static PetscTruth called = PETSC_FALSE; 1298 MPI_Comm comm = A->comm; 1299 int ierr; 1300 1301 PetscFunctionBegin; 1302 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1303 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1304 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1305 PetscFunctionReturn(0); 1306 } 1307 1308 EXTERN_C_BEGIN 1309 #undef __FUNC__ 1310 #define __FUNC__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1311 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1312 { 1313 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1314 int i,nz,n; 1315 1316 PetscFunctionBegin; 1317 nz = baij->maxnz; 1318 n = baij->n; 1319 for (i=0; i<nz; i++) { 1320 baij->j[i] = indices[i]; 1321 } 1322 baij->nz = nz; 1323 for (i=0; i<n; i++) { 1324 baij->ilen[i] = baij->imax[i]; 1325 } 1326 1327 PetscFunctionReturn(0); 1328 } 1329 EXTERN_C_END 1330 1331 #undef __FUNC__ 1332 #define __FUNC__ "MatSeqSBAIJSetColumnIndices" 1333 /*@ 1334 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1335 in the matrix. 1336 1337 Input Parameters: 1338 + mat - the SeqBAIJ matrix 1339 - indices - the column indices 1340 1341 Level: advanced 1342 1343 Notes: 1344 This can be called if you have precomputed the nonzero structure of the 1345 matrix and want to provide it to the matrix object to improve the performance 1346 of the MatSetValues() operation. 1347 1348 You MUST have set the correct numbers of nonzeros per row in the call to 1349 MatCreateSeqBAIJ(). 1350 1351 MUST be called before any calls to MatSetValues(); 1352 1353 @*/ 1354 int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices) 1355 { 1356 int ierr,(*f)(Mat,int *); 1357 1358 PetscFunctionBegin; 1359 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1360 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1361 if (f) { 1362 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1363 } else { 1364 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1365 } 1366 PetscFunctionReturn(0); 1367 } 1368 1369 /* -------------------------------------------------------------------*/ 1370 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1371 MatGetRow_SeqSBAIJ, 1372 MatRestoreRow_SeqSBAIJ, 1373 MatMult_SeqSBAIJ_N, 1374 MatMultAdd_SeqSBAIJ_N, 1375 MatMultTranspose_SeqSBAIJ, 1376 MatMultTransposeAdd_SeqSBAIJ, 1377 MatSolve_SeqSBAIJ_N, 1378 0, 1379 0, 1380 0, 1381 0, 1382 MatCholeskyFactor_SeqSBAIJ, 1383 0, 1384 MatTranspose_SeqSBAIJ, 1385 MatGetInfo_SeqSBAIJ, 1386 MatEqual_SeqSBAIJ, 1387 MatGetDiagonal_SeqSBAIJ, 1388 MatDiagonalScale_SeqSBAIJ, 1389 MatNorm_SeqSBAIJ, 1390 0, 1391 MatAssemblyEnd_SeqSBAIJ, 1392 0, 1393 MatSetOption_SeqSBAIJ, 1394 MatZeroEntries_SeqSBAIJ, 1395 MatZeroRows_SeqSBAIJ, 1396 0, 1397 0, 1398 MatCholeskyFactorSymbolic_SeqSBAIJ, 1399 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1400 MatGetSize_SeqSBAIJ, 1401 MatGetSize_SeqSBAIJ, 1402 MatGetOwnershipRange_SeqSBAIJ, 1403 0, 1404 MatIncompleteCholeskyFactorSymbolic_SeqSBAIJ, 1405 0, 1406 0, 1407 MatDuplicate_SeqSBAIJ, 1408 0, 1409 0, 1410 0, 1411 0, 1412 0, 1413 MatGetSubMatrices_SeqSBAIJ, 1414 MatIncreaseOverlap_SeqSBAIJ, 1415 MatGetValues_SeqSBAIJ, 1416 0, 1417 MatPrintHelp_SeqSBAIJ, 1418 MatScale_SeqSBAIJ, 1419 0, 1420 0, 1421 0, 1422 MatGetBlockSize_SeqSBAIJ, 1423 MatGetRowIJ_SeqSBAIJ, 1424 MatRestoreRowIJ_SeqSBAIJ, 1425 0, 1426 0, 1427 0, 1428 0, 1429 0, 1430 0, 1431 MatSetValuesBlocked_SeqSBAIJ, 1432 MatGetSubMatrix_SeqSBAIJ, 1433 0, 1434 0, 1435 MatGetMaps_Petsc}; 1436 1437 EXTERN_C_BEGIN 1438 #undef __FUNC__ 1439 #define __FUNC__ "MatStoreValues_SeqSBAIJ" 1440 int MatStoreValues_SeqSBAIJ(Mat mat) 1441 { 1442 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1443 int nz = aij->i[aij->m]*aij->bs*aij->bs2; 1444 int ierr; 1445 1446 PetscFunctionBegin; 1447 if (aij->nonew != 1) { 1448 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1449 } 1450 1451 /* allocate space for values if not already there */ 1452 if (!aij->saved_values) { 1453 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1454 } 1455 1456 /* copy values over */ 1457 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1458 PetscFunctionReturn(0); 1459 } 1460 EXTERN_C_END 1461 1462 EXTERN_C_BEGIN 1463 #undef __FUNC__ 1464 #define __FUNC__ "MatRetrieveValues_SeqSBAIJ" 1465 int MatRetrieveValues_SeqSBAIJ(Mat mat) 1466 { 1467 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1468 int nz = aij->i[aij->m]*aij->bs*aij->bs2,ierr; 1469 1470 PetscFunctionBegin; 1471 if (aij->nonew != 1) { 1472 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1473 } 1474 if (!aij->saved_values) { 1475 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1476 } 1477 1478 /* copy values over */ 1479 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 EXTERN_C_END 1483 1484 #undef __FUNC__ 1485 #define __FUNC__ "MatCreateSeqSBAIJ" 1486 /*@C 1487 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1488 compressed row) format. For good matrix assembly performance the 1489 user should preallocate the matrix storage by setting the parameter nz 1490 (or the array nnz). By setting these parameters accurately, performance 1491 during matrix assembly can be increased by more than a factor of 50. 1492 1493 Collective on MPI_Comm 1494 1495 Input Parameters: 1496 + comm - MPI communicator, set to PETSC_COMM_SELF 1497 . bs - size of block 1498 . m - number of rows, or number of columns 1499 . nz - number of block nonzeros per block row (same for all rows) 1500 - nnz - array containing the number of block nonzeros in the various block rows 1501 (possibly different for each block row) or PETSC_NULL 1502 1503 Output Parameter: 1504 . A - the symmetric matrix 1505 1506 Options Database Keys: 1507 . -mat_no_unroll - uses code that does not unroll the loops in the 1508 block calculations (much slower) 1509 . -mat_block_size - size of the blocks to use 1510 1511 Level: intermediate 1512 1513 Notes: 1514 The block AIJ format is fully compatible with standard Fortran 77 1515 storage. That is, the stored row and column indices can begin at 1516 either one (as in Fortran) or zero. See the users' manual for details. 1517 1518 Specify the preallocated storage with either nz or nnz (not both). 1519 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1520 allocation. For additional details, see the users manual chapter on 1521 matrices. 1522 1523 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1524 @*/ 1525 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1526 { 1527 Mat B; 1528 Mat_SeqSBAIJ *b; 1529 int i,len,ierr,mbs,bs2,size; 1530 PetscTruth flg; 1531 int s_nz; 1532 1533 PetscFunctionBegin; 1534 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1535 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1536 1537 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1538 mbs = m/bs; 1539 bs2 = bs*bs; 1540 1541 if (mbs*bs!=m) { 1542 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1543 } 1544 1545 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"nz cannot be less than -2: value %d",nz); 1546 if (nnz) { 1547 for (i=0; i<mbs; i++) { 1548 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1549 } 1550 } 1551 1552 *A = 0; 1553 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQSBAIJ,"Mat",comm,MatDestroy,MatView); 1554 PLogObjectCreate(B); 1555 B->data = (void*)(b = PetscNew(Mat_SeqSBAIJ));CHKPTRQ(b); 1556 ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));CHKERRQ(ierr); 1557 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1558 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1559 if (!flg) { 1560 switch (bs) { 1561 case 1: 1562 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1563 B->ops->solve = MatSolve_SeqSBAIJ_1; 1564 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1565 B->ops->mult = MatMult_SeqSBAIJ_1; 1566 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1567 break; 1568 case 2: 1569 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1570 B->ops->solve = MatSolve_SeqSBAIJ_2; 1571 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1572 B->ops->mult = MatMult_SeqSBAIJ_2; 1573 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1574 break; 1575 case 3: 1576 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1577 B->ops->solve = MatSolve_SeqSBAIJ_3; 1578 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1579 B->ops->mult = MatMult_SeqSBAIJ_3; 1580 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1581 break; 1582 case 4: 1583 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1584 B->ops->solve = MatSolve_SeqSBAIJ_4; 1585 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1586 B->ops->mult = MatMult_SeqSBAIJ_4; 1587 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1588 break; 1589 case 5: 1590 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1591 B->ops->solve = MatSolve_SeqSBAIJ_5; 1592 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1593 B->ops->mult = MatMult_SeqSBAIJ_5; 1594 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1595 break; 1596 case 6: 1597 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1598 B->ops->solve = MatSolve_SeqSBAIJ_6; 1599 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1600 B->ops->mult = MatMult_SeqSBAIJ_6; 1601 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1602 break; 1603 case 7: 1604 B->ops->mult = MatMult_SeqSBAIJ_7; 1605 B->ops->solve = MatSolve_SeqSBAIJ_7; 1606 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1607 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1608 break; 1609 } 1610 } 1611 B->ops->destroy = MatDestroy_SeqSBAIJ; 1612 B->ops->view = MatView_SeqSBAIJ; 1613 B->factor = 0; 1614 B->lupivotthreshold = 1.0; 1615 B->mapping = 0; 1616 b->row = 0; 1617 b->col = 0; 1618 b->icol = 0; 1619 b->reallocs = 0; 1620 b->saved_values = 0; 1621 1622 b->m = m; 1623 B->m = m; B->M = m; 1624 b->n = m; 1625 B->n = m; B->N = m; 1626 1627 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1628 ierr = MapCreateMPI(comm,m,m,&B->cmap);CHKERRQ(ierr); 1629 1630 b->mbs = mbs; 1631 b->nbs = mbs; 1632 b->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(b->imax); 1633 if (!nnz) { 1634 if (nz == PETSC_DEFAULT) nz = 5; 1635 else if (nz <= 0) nz = 1; 1636 for (i=0; i<mbs; i++) { 1637 b->imax[i] = (nz+1)/2; 1638 } 1639 nz = nz*mbs; 1640 } else { 1641 nz = 0; 1642 for (i=0; i<mbs; i++) {b->imax[i] = (nnz[i]+1)/2; nz += nnz[i];} 1643 } 1644 s_nz=(nz+mbs)/2; /* total diagonal and superdiagonal nonzero blocks */ 1645 1646 /* allocate the matrix space */ 1647 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (b->m+1)*sizeof(int); 1648 b->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(b->a); 1649 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1650 b->j = (int*)(b->a + s_nz*bs2); 1651 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1652 b->i = b->j + s_nz; 1653 b->singlemalloc = PETSC_TRUE; 1654 1655 /* pointer to beginning of each row */ 1656 b->i[0] = 0; 1657 for (i=1; i<mbs+1; i++) { 1658 b->i[i] = b->i[i-1] + b->imax[i-1]; 1659 } 1660 1661 1662 /* b->ilen will count nonzeros in each block row so far. */ 1663 b->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int)); 1664 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1665 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1666 1667 b->bs = bs; 1668 b->bs2 = bs2; 1669 b->mbs = mbs; 1670 b->s_nz = 0; 1671 b->s_maxnz = s_nz*bs2; 1672 b->sorted = PETSC_FALSE; 1673 b->roworiented = PETSC_TRUE; 1674 b->nonew = 0; 1675 b->diag = 0; 1676 b->solve_work = 0; 1677 b->mult_work = 0; 1678 b->spptr = 0; 1679 B->info.nz_unneeded = (PetscReal)b->s_maxnz; 1680 b->keepzeroedrows = PETSC_FALSE; 1681 1682 *A = B; 1683 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 1684 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr); } 1685 1686 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1687 "MatStoreValues_SeqSBAIJ", 1688 (void*)MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1689 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1690 "MatRetrieveValues_SeqSBAIJ", 1691 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1692 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1693 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1694 (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1695 /* printf("In MatCreateSeqSBAIJ, \n"); 1696 for (i=0; i<mbs; i++){ 1697 printf("imax[%d]= %d, ilen[%d]= %d\n", i,b->imax[i], i,b->ilen[i]); 1698 } */ 1699 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNC__ 1704 #define __FUNC__ "MatDuplicate_SeqSBAIJ" 1705 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1706 { 1707 Mat C; 1708 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1709 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1710 1711 PetscFunctionBegin; 1712 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1713 1714 *B = 0; 1715 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQSBAIJ,"Mat",A->comm,MatDestroy,MatView); 1716 PLogObjectCreate(C); 1717 C->data = (void*)(c = PetscNew(Mat_SeqSBAIJ));CHKPTRQ(c); 1718 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1719 C->ops->destroy = MatDestroy_SeqSBAIJ; 1720 C->ops->view = MatView_SeqSBAIJ; 1721 C->factor = A->factor; 1722 c->row = 0; 1723 c->col = 0; 1724 c->icol = 0; 1725 c->saved_values = 0; 1726 c->keepzeroedrows = a->keepzeroedrows; 1727 C->assembled = PETSC_TRUE; 1728 1729 c->m = C->m = a->m; 1730 c->n = C->n = a->n; 1731 C->M = a->m; 1732 C->N = a->n; 1733 1734 c->bs = a->bs; 1735 c->bs2 = a->bs2; 1736 c->mbs = a->mbs; 1737 c->nbs = a->nbs; 1738 1739 c->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->imax); 1740 c->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->ilen); 1741 for (i=0; i<mbs; i++) { 1742 c->imax[i] = a->imax[i]; 1743 c->ilen[i] = a->ilen[i]; 1744 } 1745 1746 /* allocate the matrix space */ 1747 c->singlemalloc = PETSC_TRUE; 1748 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1749 c->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(c->a); 1750 c->j = (int*)(c->a + nz*bs2); 1751 c->i = c->j + nz; 1752 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1753 if (mbs > 0) { 1754 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1755 if (cpvalues == MAT_COPY_VALUES) { 1756 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1757 } else { 1758 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1759 } 1760 } 1761 1762 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1763 c->sorted = a->sorted; 1764 c->roworiented = a->roworiented; 1765 c->nonew = a->nonew; 1766 1767 if (a->diag) { 1768 c->diag = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->diag); 1769 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1770 for (i=0; i<mbs; i++) { 1771 c->diag[i] = a->diag[i]; 1772 } 1773 } else c->diag = 0; 1774 c->s_nz = a->s_nz; 1775 c->s_maxnz = a->s_maxnz; 1776 c->solve_work = 0; 1777 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1778 c->mult_work = 0; 1779 *B = C; 1780 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1781 PetscFunctionReturn(0); 1782 } 1783 1784 #undef __FUNC__ 1785 #define __FUNC__ "MatLoad_SeqSBAIJ" 1786 int MatLoad_SeqSBAIJ(Viewer viewer,MatType type,Mat *A) 1787 { 1788 Mat_SeqSBAIJ *a; 1789 Mat B; 1790 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1791 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1792 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1793 int *masked,nmask,tmp,bs2,ishift; 1794 Scalar *aa; 1795 MPI_Comm comm = ((PetscObject)viewer)->comm; 1796 1797 PetscFunctionBegin; 1798 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1799 bs2 = bs*bs; 1800 1801 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1802 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1803 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1804 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1805 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1806 M = header[1]; N = header[2]; nz = header[3]; 1807 1808 if (header[3] < 0) { 1809 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1810 } 1811 1812 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1813 1814 /* 1815 This code adds extra rows to make sure the number of rows is 1816 divisible by the blocksize 1817 */ 1818 mbs = M/bs; 1819 extra_rows = bs - M + bs*(mbs); 1820 if (extra_rows == bs) extra_rows = 0; 1821 else mbs++; 1822 if (extra_rows) { 1823 PLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1824 } 1825 1826 /* read in row lengths */ 1827 rowlengths = (int*)PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1828 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1829 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1830 1831 /* read in column indices */ 1832 jj = (int*)PetscMalloc((nz+extra_rows)*sizeof(int));CHKPTRQ(jj); 1833 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1834 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1835 1836 /* loop over row lengths determining block row lengths */ 1837 browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1838 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1839 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1840 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1841 masked = mask + mbs; 1842 rowcount = 0; nzcount = 0; 1843 for (i=0; i<mbs; i++) { 1844 nmask = 0; 1845 for (j=0; j<bs; j++) { 1846 kmax = rowlengths[rowcount]; 1847 for (k=0; k<kmax; k++) { 1848 tmp = jj[nzcount++]/bs; 1849 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1850 } 1851 rowcount++; 1852 } 1853 browlengths[i] += nmask; 1854 /* zero out the mask elements we set */ 1855 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1856 } 1857 1858 /* create our matrix */ 1859 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1860 B = *A; 1861 a = (Mat_SeqSBAIJ*)B->data; 1862 1863 /* set matrix "i" values */ 1864 a->i[0] = 0; 1865 for (i=1; i<= mbs; i++) { 1866 a->i[i] = a->i[i-1] + browlengths[i-1]; 1867 a->ilen[i-1] = browlengths[i-1]; 1868 } 1869 a->s_nz = 0; 1870 for (i=0; i<mbs; i++) a->s_nz += browlengths[i]; 1871 1872 /* read in nonzero values */ 1873 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1874 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1875 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1876 1877 /* set "a" and "j" values into matrix */ 1878 nzcount = 0; jcount = 0; 1879 for (i=0; i<mbs; i++) { 1880 nzcountb = nzcount; 1881 nmask = 0; 1882 for (j=0; j<bs; j++) { 1883 kmax = rowlengths[i*bs+j]; 1884 for (k=0; k<kmax; k++) { 1885 tmp = jj[nzcount++]/bs; 1886 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1887 } 1888 rowcount++; 1889 } 1890 /* sort the masked values */ 1891 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1892 1893 /* set "j" values into matrix */ 1894 maskcount = 1; 1895 for (j=0; j<nmask; j++) { 1896 a->j[jcount++] = masked[j]; 1897 mask[masked[j]] = maskcount++; 1898 } 1899 /* set "a" values into matrix */ 1900 ishift = bs2*a->i[i]; 1901 for (j=0; j<bs; j++) { 1902 kmax = rowlengths[i*bs+j]; 1903 for (k=0; k<kmax; k++) { 1904 tmp = jj[nzcountb]/bs ; 1905 block = mask[tmp] - 1; 1906 point = jj[nzcountb] - bs*tmp; 1907 idx = ishift + bs2*block + j + bs*point; 1908 a->a[idx] = aa[nzcountb++]; 1909 } 1910 } 1911 /* zero out the mask elements we set */ 1912 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1913 } 1914 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1915 1916 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1917 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1918 ierr = PetscFree(aa);CHKERRQ(ierr); 1919 ierr = PetscFree(jj);CHKERRQ(ierr); 1920 ierr = PetscFree(mask);CHKERRQ(ierr); 1921 1922 B->assembled = PETSC_TRUE; 1923 1924 ierr = MatView_Private(B);CHKERRQ(ierr); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 1929 1930