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