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