1 /*$Id: sbaij.c,v 1.17 2000/08/28 16:36:39 hzhang Exp hzhang $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "petscsys.h" /*I "petscmat.h" I*/ 8 #include "src/mat/impls/baij/seq/baij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 #include "src/mat/impls/sbaij/seq/sbaij.h" 12 13 #define CHUNKSIZE 10 14 15 /* 16 Checks for missing diagonals 17 */ 18 #undef __FUNC__ 19 #define __FUNC__ "MatMissingDiagonal_SeqSBAIJ" 20 int MatMissingDiagonal_SeqSBAIJ(Mat A) 21 { 22 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 23 int *diag,*jj = a->j,i,ierr; 24 25 PetscFunctionBegin; 26 ierr = MatMarkDiagonal_SeqBAIJ(A);CHKERRQ(ierr); 27 diag = a->diag; 28 for (i=0; i<a->mbs; i++) { 29 if (jj[diag[i]] != i) { 30 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 31 } 32 } 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNC__ 37 #define __FUNC__ "MatMarkDiagonal_SeqSBAIJ" 38 int MatMarkDiagonal_SeqSBAIJ(Mat A) 39 { 40 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 41 int i,j,*diag,m = a->mbs; 42 43 PetscFunctionBegin; 44 if (a->diag) PetscFunctionReturn(0); 45 46 diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(diag); 47 PLogObjectMemory(A,(m+1)*sizeof(int)); 48 for (i=0; i<m; i++) { 49 diag[i] = a->i[i+1]; 50 for (j=a->i[i]; j<a->i[i+1]; j++) { 51 if (a->j[j] == i) { 52 diag[i] = j; 53 break; 54 } 55 } 56 } 57 a->diag = diag; 58 PetscFunctionReturn(0); 59 } 60 61 62 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 63 64 #undef __FUNC__ 65 #define __FUNC__ "MatGetRowIJ_SeqSBAIJ" 66 static int MatGetRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 67 { 68 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 69 int ierr,n = a->mbs,i; 70 71 PetscFunctionBegin; 72 *nn = n; 73 if (!ia) PetscFunctionReturn(0); 74 if (symmetric) { 75 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 76 } else if (oshift == 1) { 77 /* temporarily add 1 to i and j indices */ 78 int nz = a->i[n] + 1; 79 for (i=0; i<nz; i++) a->j[i]++; 80 for (i=0; i<n+1; i++) a->i[i]++; 81 *ia = a->i; *ja = a->j; 82 } else { 83 *ia = a->i; *ja = a->j; 84 } 85 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNC__ 90 #define __FUNC__ "MatRestoreRowIJ_SeqSBAIJ" 91 static int MatRestoreRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 92 PetscTruth *done) 93 { 94 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 95 int i,n = a->mbs,ierr; 96 97 PetscFunctionBegin; 98 if (!ia) PetscFunctionReturn(0); 99 if (symmetric) { 100 ierr = PetscFree(*ia);CHKERRQ(ierr); 101 ierr = PetscFree(*ja);CHKERRQ(ierr); 102 } else if (oshift == 1) { 103 int nz = a->i[n]; 104 for (i=0; i<nz; i++) a->j[i]--; 105 for (i=0; i<n+1; i++) a->i[i]--; 106 } 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNC__ 111 #define __FUNC__ "MatGetBlockSize_SeqSBAIJ" 112 int MatGetBlockSize_SeqSBAIJ(Mat mat,int *bs) 113 { 114 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)mat->data; 115 116 PetscFunctionBegin; 117 *bs = sbaij->bs; 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNC__ 122 #define __FUNC__ "MatDestroy_SeqSBAIJ" 123 int MatDestroy_SeqSBAIJ(Mat A) 124 { 125 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 126 int ierr; 127 128 PetscFunctionBegin; 129 if (A->mapping) { 130 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 131 } 132 if (A->bmapping) { 133 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 134 } 135 if (A->rmap) { 136 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 137 } 138 if (A->cmap) { 139 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 140 } 141 #if defined(PETSC_USE_LOG) 142 PLogObjectState((PetscObject)A,"Rows=%d, s_NZ=%d",a->m,a->s_nz); 143 #endif 144 ierr = PetscFree(a->a);CHKERRQ(ierr); 145 if (!a->singlemalloc) { 146 ierr = PetscFree(a->i);CHKERRQ(ierr); 147 ierr = PetscFree(a->j);CHKERRQ(ierr); 148 } 149 if (a->row) { 150 ierr = ISDestroy(a->row);CHKERRQ(ierr); 151 } 152 if (a->col) { 153 ierr = ISDestroy(a->col);CHKERRQ(ierr); 154 } 155 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 156 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 157 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 158 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 159 if (a->mult_work) {ierr = PetscFree(a->mult_work);CHKERRQ(ierr);} 160 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 161 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 162 ierr = PetscFree(a);CHKERRQ(ierr); 163 PLogObjectDestroy(A); 164 PetscHeaderDestroy(A); 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNC__ 169 #define __FUNC__ "MatSetOption_SeqSBAIJ" 170 int MatSetOption_SeqSBAIJ(Mat A,MatOption op) 171 { 172 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 173 174 PetscFunctionBegin; 175 if (op == MAT_ROW_ORIENTED) a->roworiented = PETSC_TRUE; 176 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = PETSC_FALSE; 177 else if (op == MAT_COLUMNS_SORTED) a->sorted = PETSC_TRUE; 178 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = PETSC_FALSE; 179 else if (op == MAT_KEEP_ZEROED_ROWS) a->keepzeroedrows = PETSC_TRUE; 180 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 181 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 182 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 183 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 184 else if (op == MAT_ROWS_SORTED || 185 op == MAT_ROWS_UNSORTED || 186 op == MAT_SYMMETRIC || 187 op == MAT_STRUCTURALLY_SYMMETRIC || 188 op == MAT_YES_NEW_DIAGONALS || 189 op == MAT_IGNORE_OFF_PROC_ENTRIES || 190 op == MAT_USE_HASH_TABLE) { 191 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 192 } else if (op == MAT_NO_NEW_DIAGONALS) { 193 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 194 } else { 195 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 196 } 197 PetscFunctionReturn(0); 198 } 199 200 201 #undef __FUNC__ 202 #define __FUNC__ "MatGetSize_SeqSBAIJ" 203 int MatGetSize_SeqSBAIJ(Mat A,int *m,int *n) 204 { 205 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 206 207 PetscFunctionBegin; 208 if (m) *m = a->m; 209 if (n) *n = a->m; 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNC__ 214 #define __FUNC__ "MatGetOwnershipRange_SeqSBAIJ" 215 int MatGetOwnershipRange_SeqSBAIJ(Mat A,int *m,int *n) 216 { 217 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 218 219 PetscFunctionBegin; 220 *m = 0; *n = a->m; 221 PetscFunctionReturn(0); 222 } 223 224 #undef __FUNC__ 225 #define __FUNC__ "MatGetRow_SeqSBAIJ" 226 int MatGetRow_SeqSBAIJ(Mat A,int row,int *ncols,int **cols,Scalar **v) 227 { 228 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 229 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 230 MatScalar *aa,*aa_i; 231 Scalar *v_i; 232 233 PetscFunctionBegin; 234 bs = a->bs; 235 ai = a->i; 236 aj = a->j; 237 aa = a->a; 238 bs2 = a->bs2; 239 240 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 241 242 bn = row/bs; /* Block number */ 243 bp = row % bs; /* Block position */ 244 M = ai[bn+1] - ai[bn]; 245 *ncols = bs*M; 246 247 if (v) { 248 *v = 0; 249 if (*ncols) { 250 *v = (Scalar*)PetscMalloc((*ncols+row)*sizeof(Scalar));CHKPTRQ(*v); 251 for (i=0; i<M; i++) { /* for each block in the block row */ 252 v_i = *v + i*bs; 253 aa_i = aa + bs2*(ai[bn] + i); 254 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 255 } 256 } 257 } 258 259 if (cols) { 260 *cols = 0; 261 if (*ncols) { 262 *cols = (int*)PetscMalloc((*ncols+row)*sizeof(int));CHKPTRQ(*cols); 263 for (i=0; i<M; i++) { /* for each block in the block row */ 264 cols_i = *cols + i*bs; 265 itmp = bs*aj[ai[bn] + i]; 266 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 267 } 268 } 269 } 270 271 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 272 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 273 #ifdef column_search 274 v_i = *v + M*bs; 275 cols_i = *cols + M*bs; 276 for (i=0; i<bn; i++){ /* for each block row */ 277 M = ai[i+1] - ai[i]; 278 for (j=0; j<M; j++){ 279 itmp = aj[ai[i] + j]; /* block column value */ 280 if (itmp == bn){ 281 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 282 for (k=0; k<bs; k++) { 283 *cols_i++ = i*bs+k; 284 *v_i++ = aa_i[k]; 285 } 286 *ncols += bs; 287 break; 288 } 289 } 290 } 291 #endif 292 293 PetscFunctionReturn(0); 294 } 295 296 #undef __FUNC__ 297 #define __FUNC__ "MatRestoreRow_SeqSBAIJ" 298 int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 299 { 300 int ierr; 301 302 PetscFunctionBegin; 303 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 304 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 305 PetscFunctionReturn(0); 306 } 307 308 #undef __FUNC__ 309 #define __FUNC__ "MatTranspose_SeqSBAIJ" 310 int MatTranspose_SeqSBAIJ(Mat A,Mat *B) 311 { 312 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 313 Mat C; 314 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 315 int *rows,*cols,bs2=a->bs2,refct; 316 MatScalar *array = a->a; 317 318 PetscFunctionBegin; 319 if (!B && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 320 col = (int*)PetscMalloc((1+nbs)*sizeof(int));CHKPTRQ(col); 321 ierr = PetscMemzero(col,(1+nbs)*sizeof(int));CHKERRQ(ierr); 322 323 for (i=0; i<ai[mbs]; i++) col[aj[i]] += 1; 324 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C);CHKERRQ(ierr); 325 ierr = PetscFree(col);CHKERRQ(ierr); 326 rows = (int*)PetscMalloc(2*bs*sizeof(int));CHKPTRQ(rows); 327 cols = rows + bs; 328 for (i=0; i<mbs; i++) { 329 cols[0] = i*bs; 330 for (k=1; k<bs; k++) cols[k] = cols[k-1] + 1; 331 len = ai[i+1] - ai[i]; 332 for (j=0; j<len; j++) { 333 rows[0] = (*aj++)*bs; 334 for (k=1; k<bs; k++) rows[k] = rows[k-1] + 1; 335 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES);CHKERRQ(ierr); 336 array += bs2; 337 } 338 } 339 ierr = PetscFree(rows);CHKERRQ(ierr); 340 341 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 342 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 343 344 if (B) { 345 *B = C; 346 } else { 347 PetscOps *Abops; 348 MatOps Aops; 349 350 /* This isn't really an in-place transpose */ 351 ierr = PetscFree(a->a);CHKERRQ(ierr); 352 if (!a->singlemalloc) { 353 ierr = PetscFree(a->i);CHKERRQ(ierr); 354 ierr = PetscFree(a->j);CHKERRQ(ierr); 355 } 356 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 357 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 358 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 359 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 360 ierr = PetscFree(a);CHKERRQ(ierr); 361 362 363 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 364 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 365 366 /* 367 This is horrible, horrible code. We need to keep the 368 A pointers for the bops and ops but copy everything 369 else from C. 370 */ 371 Abops = A->bops; 372 Aops = A->ops; 373 refct = A->refct; 374 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 375 A->bops = Abops; 376 A->ops = Aops; 377 A->qlist = 0; 378 A->refct = refct; 379 380 PetscHeaderDestroy(C); 381 } 382 PetscFunctionReturn(0); 383 } 384 385 #undef __FUNC__ 386 #define __FUNC__ "MatView_SeqSBAIJ_Binary" 387 static int MatView_SeqSBAIJ_Binary(Mat A,Viewer viewer) 388 { 389 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 390 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 391 Scalar *aa; 392 FILE *file; 393 394 PetscFunctionBegin; 395 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 396 col_lens = (int*)PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 397 col_lens[0] = MAT_COOKIE; 398 399 col_lens[1] = a->m; 400 col_lens[2] = a->m; 401 col_lens[3] = a->s_nz*bs2; 402 403 /* store lengths of each row and write (including header) to file */ 404 count = 0; 405 for (i=0; i<a->mbs; i++) { 406 for (j=0; j<bs; j++) { 407 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 408 } 409 } 410 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1);CHKERRQ(ierr); 411 ierr = PetscFree(col_lens);CHKERRQ(ierr); 412 413 /* store column indices (zero start index) */ 414 jj = (int*)PetscMalloc((a->s_nz+1)*bs2*sizeof(int));CHKPTRQ(jj); 415 count = 0; 416 for (i=0; i<a->mbs; i++) { 417 for (j=0; j<bs; j++) { 418 for (k=a->i[i]; k<a->i[i+1]; k++) { 419 for (l=0; l<bs; l++) { 420 jj[count++] = bs*a->j[k] + l; 421 } 422 } 423 } 424 } 425 ierr = PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);CHKERRQ(ierr); 426 ierr = PetscFree(jj);CHKERRQ(ierr); 427 428 /* store nonzero values */ 429 aa = (Scalar*)PetscMalloc((a->s_nz+1)*bs2*sizeof(Scalar));CHKPTRQ(aa); 430 count = 0; 431 for (i=0; i<a->mbs; i++) { 432 for (j=0; j<bs; j++) { 433 for (k=a->i[i]; k<a->i[i+1]; k++) { 434 for (l=0; l<bs; l++) { 435 aa[count++] = a->a[bs2*k + l*bs + j]; 436 } 437 } 438 } 439 } 440 ierr = PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);CHKERRQ(ierr); 441 ierr = PetscFree(aa);CHKERRQ(ierr); 442 443 ierr = ViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 444 if (file) { 445 fprintf(file,"-matload_block_size %d\n",a->bs); 446 } 447 PetscFunctionReturn(0); 448 } 449 450 #undef __FUNC__ 451 #define __FUNC__ "MatView_SeqSBAIJ_ASCII" 452 static int MatView_SeqSBAIJ_ASCII(Mat A,Viewer viewer) 453 { 454 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 455 int ierr,i,j,format,bs = a->bs,k,l,bs2=a->bs2; 456 char *outputname; 457 458 PetscFunctionBegin; 459 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 460 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 461 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 462 ierr = ViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 463 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 464 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 465 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 466 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 467 for (i=0; i<a->mbs; i++) { 468 for (j=0; j<bs; j++) { 469 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 470 for (k=a->i[i]; k<a->i[i+1]; k++) { 471 for (l=0; l<bs; l++) { 472 #if defined(PETSC_USE_COMPLEX) 473 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 474 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 475 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 476 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 477 ierr = ViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 478 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 479 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 480 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 481 } 482 #else 483 if (a->a[bs2*k + l*bs + j] != 0.0) { 484 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 485 } 486 #endif 487 } 488 } 489 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 490 } 491 } 492 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 493 } else { 494 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 495 for (i=0; i<a->mbs; i++) { 496 for (j=0; j<bs; j++) { 497 ierr = ViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 498 for (k=a->i[i]; k<a->i[i+1]; k++) { 499 for (l=0; l<bs; l++) { 500 #if defined(PETSC_USE_COMPLEX) 501 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 502 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 503 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 504 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 505 ierr = ViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 506 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 507 } else { 508 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 509 } 510 #else 511 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 512 #endif 513 } 514 } 515 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 516 } 517 } 518 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 519 } 520 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 524 #undef __FUNC__ 525 #define __FUNC__ "MatView_SeqSBAIJ_Draw_Zoom" 526 static int MatView_SeqSBAIJ_Draw_Zoom(Draw draw,void *Aa) 527 { 528 Mat A = (Mat) Aa; 529 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 530 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 531 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 532 MatScalar *aa; 533 MPI_Comm comm; 534 Viewer viewer; 535 536 PetscFunctionBegin; 537 /* 538 This is nasty. If this is called from an originally parallel matrix 539 then all processes call this,but only the first has the matrix so the 540 rest should return immediately. 541 */ 542 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 543 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 544 if (rank) PetscFunctionReturn(0); 545 546 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 547 548 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 549 DrawString(draw, .3*(xl+xr), .3*(yl+yr), DRAW_BLACK, "symmetric"); 550 551 /* loop over matrix elements drawing boxes */ 552 color = DRAW_BLUE; 553 for (i=0,row=0; i<mbs; i++,row+=bs) { 554 for (j=a->i[i]; j<a->i[i+1]; j++) { 555 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 556 x_l = a->j[j]*bs; x_r = x_l + 1.0; 557 aa = a->a + j*bs2; 558 for (k=0; k<bs; k++) { 559 for (l=0; l<bs; l++) { 560 if (PetscRealPart(*aa++) >= 0.) continue; 561 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 562 } 563 } 564 } 565 } 566 color = DRAW_CYAN; 567 for (i=0,row=0; i<mbs; i++,row+=bs) { 568 for (j=a->i[i]; j<a->i[i+1]; j++) { 569 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 570 x_l = a->j[j]*bs; x_r = x_l + 1.0; 571 aa = a->a + j*bs2; 572 for (k=0; k<bs; k++) { 573 for (l=0; l<bs; l++) { 574 if (PetscRealPart(*aa++) != 0.) continue; 575 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 576 } 577 } 578 } 579 } 580 581 color = DRAW_RED; 582 for (i=0,row=0; i<mbs; i++,row+=bs) { 583 for (j=a->i[i]; j<a->i[i+1]; j++) { 584 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 585 x_l = a->j[j]*bs; x_r = x_l + 1.0; 586 aa = a->a + j*bs2; 587 for (k=0; k<bs; k++) { 588 for (l=0; l<bs; l++) { 589 if (PetscRealPart(*aa++) <= 0.) continue; 590 ierr = DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 591 } 592 } 593 } 594 } 595 PetscFunctionReturn(0); 596 } 597 598 #undef __FUNC__ 599 #define __FUNC__ "MatView_SeqSBAIJ_Draw" 600 static int MatView_SeqSBAIJ_Draw(Mat A,Viewer viewer) 601 { 602 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 603 int ierr; 604 PetscReal xl,yl,xr,yr,w,h; 605 Draw draw; 606 PetscTruth isnull; 607 608 PetscFunctionBegin; 609 610 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 611 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 612 613 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 614 xr = a->m; yr = a->m; h = yr/10.0; w = xr/10.0; 615 xr += w; yr += h; xl = -w; yl = -h; 616 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 617 ierr = DrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 618 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 619 PetscFunctionReturn(0); 620 } 621 622 #undef __FUNC__ 623 #define __FUNC__ "MatView_SeqSBAIJ" 624 int MatView_SeqSBAIJ(Mat A,Viewer viewer) 625 { 626 int ierr; 627 PetscTruth issocket,isascii,isbinary,isdraw; 628 629 PetscFunctionBegin; 630 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 631 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 632 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 633 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 634 if (issocket) { 635 SETERRQ(PETSC_ERR_SUP,0,"Socket viewer not supported"); 636 } else if (isascii){ 637 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 638 } else if (isbinary) { 639 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 640 } else if (isdraw) { 641 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 642 } else { 643 SETERRQ1(1,1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 644 } 645 PetscFunctionReturn(0); 646 } 647 648 649 #undef __FUNC__ 650 #define __FUNC__ "MatGetValues_SeqSBAIJ" 651 int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 652 { 653 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 654 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 655 int *ai = a->i,*ailen = a->ilen; 656 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 657 MatScalar *ap,*aa = a->a,zero = 0.0; 658 659 PetscFunctionBegin; 660 for (k=0; k<m; k++) { /* loop over rows */ 661 row = im[k]; brow = row/bs; 662 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 663 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 664 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 665 nrow = ailen[brow]; 666 for (l=0; l<n; l++) { /* loop over columns */ 667 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 668 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 669 col = in[l] ; 670 bcol = col/bs; 671 cidx = col%bs; 672 ridx = row%bs; 673 high = nrow; 674 low = 0; /* assume unsorted */ 675 while (high-low > 5) { 676 t = (low+high)/2; 677 if (rp[t] > bcol) high = t; 678 else low = t; 679 } 680 for (i=low; i<high; i++) { 681 if (rp[i] > bcol) break; 682 if (rp[i] == bcol) { 683 *v++ = ap[bs2*i+bs*cidx+ridx]; 684 goto finished; 685 } 686 } 687 *v++ = zero; 688 finished:; 689 } 690 } 691 PetscFunctionReturn(0); 692 } 693 694 695 #undef __FUNC__ 696 #define __FUNC__ "MatSetValuesBlocked_SeqSBAIJ" 697 int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 698 { 699 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 700 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 701 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 702 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval,ierr; 703 Scalar *value = v; 704 MatScalar *ap,*aa=a->a,*bap; 705 706 PetscFunctionBegin; 707 if (roworiented) { 708 stepval = (n-1)*bs; 709 } else { 710 stepval = (m-1)*bs; 711 } 712 for (k=0; k<m; k++) { /* loop over added rows */ 713 row = im[k]; 714 if (row < 0) continue; 715 #if defined(PETSC_USE_BOPT_g) 716 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 717 #endif 718 rp = aj + ai[row]; 719 ap = aa + bs2*ai[row]; 720 rmax = imax[row]; 721 nrow = ailen[row]; 722 low = 0; 723 for (l=0; l<n; l++) { /* loop over added columns */ 724 if (in[l] < 0) continue; 725 #if defined(PETSC_USE_BOPT_g) 726 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 727 #endif 728 col = in[l]; 729 if (roworiented) { 730 value = v + k*(stepval+bs)*bs + l*bs; 731 } else { 732 value = v + l*(stepval+bs)*bs + k*bs; 733 } 734 if (!sorted) low = 0; high = nrow; 735 while (high-low > 7) { 736 t = (low+high)/2; 737 if (rp[t] > col) high = t; 738 else low = t; 739 } 740 for (i=low; i<high; i++) { 741 if (rp[i] > col) break; 742 if (rp[i] == col) { 743 bap = ap + bs2*i; 744 if (roworiented) { 745 if (is == ADD_VALUES) { 746 for (ii=0; ii<bs; ii++,value+=stepval) { 747 for (jj=ii; jj<bs2; jj+=bs) { 748 bap[jj] += *value++; 749 } 750 } 751 } else { 752 for (ii=0; ii<bs; ii++,value+=stepval) { 753 for (jj=ii; jj<bs2; jj+=bs) { 754 bap[jj] = *value++; 755 } 756 } 757 } 758 } else { 759 if (is == ADD_VALUES) { 760 for (ii=0; ii<bs; ii++,value+=stepval) { 761 for (jj=0; jj<bs; jj++) { 762 *bap++ += *value++; 763 } 764 } 765 } else { 766 for (ii=0; ii<bs; ii++,value+=stepval) { 767 for (jj=0; jj<bs; jj++) { 768 *bap++ = *value++; 769 } 770 } 771 } 772 } 773 goto noinsert2; 774 } 775 } 776 if (nonew == 1) goto noinsert2; 777 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 778 if (nrow >= rmax) { 779 /* there is no extra room in row, therefore enlarge */ 780 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 781 MatScalar *new_a; 782 783 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 784 785 /* malloc new storage space */ 786 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 787 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 788 new_j = (int*)(new_a + bs2*new_nz); 789 new_i = new_j + new_nz; 790 791 /* copy over old data into new slots */ 792 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 793 for (ii=row+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 794 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 795 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 796 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 797 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 798 ierr = PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 799 ierr = PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 800 /* free up old matrix storage */ 801 ierr = PetscFree(a->a);CHKERRQ(ierr); 802 if (!a->singlemalloc) { 803 ierr = PetscFree(a->i);CHKERRQ(ierr); 804 ierr = PetscFree(a->j);CHKERRQ(ierr); 805 } 806 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 807 a->singlemalloc = PETSC_TRUE; 808 809 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 810 rmax = imax[row] = imax[row] + CHUNKSIZE; 811 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 812 a->maxnz += bs2*CHUNKSIZE; 813 a->reallocs++; 814 a->nz++; 815 } 816 N = nrow++ - 1; 817 /* shift up all the later entries in this row */ 818 for (ii=N; ii>=i; ii--) { 819 rp[ii+1] = rp[ii]; 820 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 821 } 822 if (N >= i) { 823 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 824 } 825 rp[i] = col; 826 bap = ap + bs2*i; 827 if (roworiented) { 828 for (ii=0; ii<bs; ii++,value+=stepval) { 829 for (jj=ii; jj<bs2; jj+=bs) { 830 bap[jj] = *value++; 831 } 832 } 833 } else { 834 for (ii=0; ii<bs; ii++,value+=stepval) { 835 for (jj=0; jj<bs; jj++) { 836 *bap++ = *value++; 837 } 838 } 839 } 840 noinsert2:; 841 low = i; 842 } 843 ailen[row] = nrow; 844 } 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNC__ 849 #define __FUNC__ "MatAssemblyEnd_SeqSBAIJ" 850 int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 851 { 852 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 853 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 854 int m = a->m,*ip,N,*ailen = a->ilen; 855 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 856 MatScalar *aa = a->a,*ap; 857 858 PetscFunctionBegin; 859 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 860 861 if (m) rmax = ailen[0]; 862 for (i=1; i<mbs; i++) { 863 /* move each row back by the amount of empty slots (fshift) before it*/ 864 fshift += imax[i-1] - ailen[i-1]; 865 rmax = PetscMax(rmax,ailen[i]); 866 if (fshift) { 867 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 868 N = ailen[i]; 869 for (j=0; j<N; j++) { 870 ip[j-fshift] = ip[j]; 871 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 872 } 873 } 874 ai[i] = ai[i-1] + ailen[i-1]; 875 } 876 if (mbs) { 877 fshift += imax[mbs-1] - ailen[mbs-1]; 878 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 879 } 880 /* reset ilen and imax for each row */ 881 for (i=0; i<mbs; i++) { 882 ailen[i] = imax[i] = ai[i+1] - ai[i]; 883 } 884 a->s_nz = ai[mbs]; 885 886 /* diagonals may have moved, so kill the diagonal pointers */ 887 if (fshift && a->diag) { 888 ierr = PetscFree(a->diag);CHKERRQ(ierr); 889 PLogObjectMemory(A,-(m+1)*sizeof(int)); 890 a->diag = 0; 891 } 892 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 893 m,a->m,a->bs,fshift*bs2,a->s_nz*bs2); 894 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 895 a->reallocs); 896 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 897 a->reallocs = 0; 898 A->info.nz_unneeded = (PetscReal)fshift*bs2; 899 900 PetscFunctionReturn(0); 901 } 902 903 /* 904 This function returns an array of flags which indicate the locations of contiguous 905 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 906 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 907 Assume: sizes should be long enough to hold all the values. 908 */ 909 #undef __FUNC__ 910 #define __FUNC__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 911 static int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 912 { 913 int i,j,k,row; 914 PetscTruth flg; 915 916 PetscFunctionBegin; 917 for (i=0,j=0; i<n; j++) { 918 row = idx[i]; 919 if (row%bs!=0) { /* Not the begining of a block */ 920 sizes[j] = 1; 921 i++; 922 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 923 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 924 i++; 925 } else { /* Begining of the block, so check if the complete block exists */ 926 flg = PETSC_TRUE; 927 for (k=1; k<bs; k++) { 928 if (row+k != idx[i+k]) { /* break in the block */ 929 flg = PETSC_FALSE; 930 break; 931 } 932 } 933 if (flg == PETSC_TRUE) { /* No break in the bs */ 934 sizes[j] = bs; 935 i+= bs; 936 } else { 937 sizes[j] = 1; 938 i++; 939 } 940 } 941 } 942 *bs_max = j; 943 PetscFunctionReturn(0); 944 } 945 946 #undef __FUNC__ 947 #define __FUNC__ "MatZeroRows_SeqSBAIJ" 948 int MatZeroRows_SeqSBAIJ(Mat A,IS is,Scalar *diag) 949 { 950 Mat_SeqSBAIJ *sbaij=(Mat_SeqSBAIJ*)A->data; 951 int ierr,i,j,k,count,is_n,*is_idx,*rows; 952 int bs=sbaij->bs,bs2=sbaij->bs2,*sizes,row,bs_max; 953 Scalar zero = 0.0; 954 MatScalar *aa; 955 956 PetscFunctionBegin; 957 /* Make a copy of the IS and sort it */ 958 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 959 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 960 961 /* allocate memory for rows,sizes */ 962 rows = (int*)PetscMalloc((3*is_n+1)*sizeof(int));CHKPTRQ(rows); 963 sizes = rows + is_n; 964 965 /* initialize copy IS values to rows, and sort them */ 966 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 967 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 968 if (sbaij->keepzeroedrows) { /* do not change nonzero structure */ 969 for (i=0; i<is_n; i++) { sizes[i] = 1; } /* sizes: size of blocks, = 1 or bs */ 970 bs_max = is_n; /* bs_max: num. of contiguous block row in the row */ 971 } else { 972 ierr = MatZeroRows_SeqSBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 973 } 974 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 975 976 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 977 row = rows[j]; /* row to be zeroed */ 978 if (row < 0 || row > sbaij->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"row %d out of range",row); 979 count = (sbaij->i[row/bs +1] - sbaij->i[row/bs])*bs; /* num. of elements in the row */ 980 aa = sbaij->a + sbaij->i[row/bs]*bs2 + (row%bs); 981 if (sizes[i] == bs && !sbaij->keepzeroedrows) { 982 if (diag) { 983 if (sbaij->ilen[row/bs] > 0) { 984 sbaij->ilen[row/bs] = 1; 985 sbaij->j[sbaij->i[row/bs]] = row/bs; 986 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 987 } 988 /* Now insert all the diagoanl values for this bs */ 989 for (k=0; k<bs; k++) { 990 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 991 } 992 } else { /* (!diag) */ 993 sbaij->ilen[row/bs] = 0; 994 } /* end (!diag) */ 995 } else { /* (sizes[i] != bs), broken block */ 996 #if defined (PETSC_USE_DEBUG) 997 if (sizes[i] != 1) SETERRQ(1,0,"Internal Error. Value should be 1"); 998 #endif 999 for (k=0; k<count; k++) { 1000 aa[0] = zero; 1001 aa+=bs; 1002 } 1003 if (diag) { 1004 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 1005 } 1006 } 1007 } 1008 1009 ierr = PetscFree(rows);CHKERRQ(ierr); 1010 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1012 PetscFunctionReturn(0); 1013 } 1014 1015 /* Only add/insert a(i,j) with i<=j (blocks). 1016 Any a(i,j) with i>j input by user is ingored. 1017 */ 1018 1019 #undef __FUNC__ 1020 #define __FUNC__ "MatSetValues_SeqSBAIJ" 1021 int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 1022 { 1023 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1024 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 1025 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 1026 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 1027 int ridx,cidx,bs2=a->bs2,ierr; 1028 MatScalar *ap,value,*aa=a->a,*bap; 1029 1030 PetscFunctionBegin; 1031 1032 for (k=0; k<m; k++) { /* loop over added rows */ 1033 row = im[k]; /* row number */ 1034 brow = row/bs; /* block row number */ 1035 if (row < 0) continue; 1036 #if defined(PETSC_USE_BOPT_g) 1037 if (row >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large: row %d max %d",row,a->m); 1038 #endif 1039 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 1040 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 1041 rmax = imax[brow]; /* maximum space allocated for this row */ 1042 nrow = ailen[brow]; /* actual length of this row */ 1043 low = 0; 1044 1045 for (l=0; l<n; l++) { /* loop over added columns */ 1046 if (in[l] < 0) continue; 1047 #if defined(PETSC_USE_BOPT_g) 1048 if (in[l] >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large: col %d max %d",in[l],a->m); 1049 #endif 1050 col = in[l]; 1051 bcol = col/bs; /* block col number */ 1052 1053 if (brow <= bcol){ 1054 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 1055 /* if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ */ 1056 /* element value a(k,l) */ 1057 if (roworiented) { 1058 value = v[l + k*n]; 1059 } else { 1060 value = v[k + l*m]; 1061 } 1062 1063 /* move pointer bap to a(k,l) quickly and add/insert value */ 1064 if (!sorted) low = 0; high = nrow; 1065 while (high-low > 7) { 1066 t = (low+high)/2; 1067 if (rp[t] > bcol) high = t; 1068 else low = t; 1069 } 1070 for (i=low; i<high; i++) { 1071 /* printf("The loop of i=low.., rp[%d]=%d\n",i,rp[i]); */ 1072 if (rp[i] > bcol) break; 1073 if (rp[i] == bcol) { 1074 bap = ap + bs2*i + bs*cidx + ridx; 1075 if (is == ADD_VALUES) *bap += value; 1076 else *bap = value; 1077 goto noinsert1; 1078 } 1079 } 1080 1081 if (nonew == 1) goto noinsert1; 1082 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1083 if (nrow >= rmax) { 1084 /* there is no extra room in row, therefore enlarge */ 1085 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 1086 MatScalar *new_a; 1087 1088 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 1089 1090 /* Malloc new storage space */ 1091 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 1092 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 1093 new_j = (int*)(new_a + bs2*new_nz); 1094 new_i = new_j + new_nz; 1095 1096 /* copy over old data into new slots */ 1097 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 1098 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 1099 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 1100 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 1101 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 1102 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1103 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 1104 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 1105 /* free up old matrix storage */ 1106 ierr = PetscFree(a->a);CHKERRQ(ierr); 1107 if (!a->singlemalloc) { 1108 ierr = PetscFree(a->i);CHKERRQ(ierr); 1109 ierr = PetscFree(a->j);CHKERRQ(ierr); 1110 } 1111 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 1112 a->singlemalloc = PETSC_TRUE; 1113 1114 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1115 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1116 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 1117 a->s_maxnz += bs2*CHUNKSIZE; 1118 a->reallocs++; 1119 a->s_nz++; 1120 } 1121 1122 N = nrow++ - 1; 1123 /* shift up all the later entries in this row */ 1124 for (ii=N; ii>=i; ii--) { 1125 rp[ii+1] = rp[ii]; 1126 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1127 } 1128 if (N>=i) { 1129 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1130 } 1131 rp[i] = bcol; 1132 ap[bs2*i + bs*cidx + ridx] = value; 1133 noinsert1:; 1134 low = i; 1135 /* } */ 1136 } /* end of if .. if.. */ 1137 } /* end of loop over added columns */ 1138 ailen[brow] = nrow; 1139 } /* end of loop over added rows */ 1140 1141 PetscFunctionReturn(0); 1142 } 1143 1144 extern int MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*); 1145 extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal); 1146 extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int); 1147 extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 1148 extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 1149 extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec); 1150 extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec); 1151 extern int MatScale_SeqSBAIJ(Scalar*,Mat); 1152 extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 1153 extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 1154 extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec); 1155 extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 1156 extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 1157 extern int MatZeroEntries_SeqSBAIJ(Mat); 1158 1159 extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 1160 extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 1161 extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 1162 extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 1163 extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 1164 extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 1165 extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 1166 extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 1167 extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec); 1168 extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec); 1169 extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec); 1170 extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec); 1171 extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec); 1172 extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec); 1173 extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec); 1174 1175 extern int MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 1176 extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 1177 extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 1178 extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 1179 extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 1180 extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 1181 extern int MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,Mat*); 1182 1183 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 1184 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 1185 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 1186 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 1187 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 1188 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 1189 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 1190 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 1191 1192 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 1193 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 1194 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 1195 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 1196 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 1197 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 1198 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 1199 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 1200 1201 #ifdef MatIncompleteCholeskyFactor 1202 /* This function is modified from MatILUFactor_SeqSBAIJ. 1203 Needs further work! Don't forget to add the function to the matrix table. */ 1204 #undef __FUNC__ 1205 #define __FUNC__ "MatIncompleteCholeskyFactor_SeqSBAIJ" 1206 int MatIncompleteCholeskyFactor_SeqSBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1207 { 1208 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1209 Mat outA; 1210 int ierr; 1211 PetscTruth row_identity,col_identity; 1212 1213 PetscFunctionBegin; 1214 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,0,"Only levels = 0 supported for in-place ILU"); 1215 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1216 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1217 if (!row_identity || !col_identity) { 1218 SETERRQ(1,1,"Row and column permutations must be identity for in-place ILU"); 1219 } 1220 1221 outA = inA; 1222 inA->factor = FACTOR_LU; 1223 1224 if (!a->diag) { 1225 ierr = MatMarkDiagonal_SeqBAIJ(inA);CHKERRQ(ierr); 1226 } 1227 /* 1228 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 1229 for ILU(0) factorization with natural ordering 1230 */ 1231 switch (a->bs) { 1232 case 1: 1233 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1234 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1235 case 2: 1236 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 1237 inA->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 1238 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1239 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1240 break; 1241 case 3: 1242 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 1243 inA->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; 1244 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; 1245 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1246 break; 1247 case 4: 1248 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1249 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1250 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1251 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1252 break; 1253 case 5: 1254 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1255 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1256 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1257 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1258 break; 1259 case 6: 1260 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1261 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1262 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1263 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1264 break; 1265 case 7: 1266 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1267 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1268 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1269 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1270 break; 1271 default: 1272 a->row = row; 1273 a->col = col; 1274 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1275 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1276 1277 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1278 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1279 PLogObjectParent(inA,a->icol); 1280 1281 if (!a->solve_work) { 1282 a->solve_work = (Scalar*)PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1283 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1284 } 1285 } 1286 1287 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1288 1289 PetscFunctionReturn(0); 1290 } 1291 #endif 1292 1293 #undef __FUNC__ 1294 #define __FUNC__ "MatPrintHelp_SeqSBAIJ" 1295 int MatPrintHelp_SeqSBAIJ(Mat A) 1296 { 1297 static PetscTruth called = PETSC_FALSE; 1298 MPI_Comm comm = A->comm; 1299 int ierr; 1300 1301 PetscFunctionBegin; 1302 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1303 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1304 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1305 PetscFunctionReturn(0); 1306 } 1307 1308 EXTERN_C_BEGIN 1309 #undef __FUNC__ 1310 #define __FUNC__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1311 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1312 { 1313 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1314 int i,nz,n; 1315 1316 PetscFunctionBegin; 1317 nz = baij->maxnz; 1318 n = baij->n; 1319 for (i=0; i<nz; i++) { 1320 baij->j[i] = indices[i]; 1321 } 1322 baij->nz = nz; 1323 for (i=0; i<n; i++) { 1324 baij->ilen[i] = baij->imax[i]; 1325 } 1326 1327 PetscFunctionReturn(0); 1328 } 1329 EXTERN_C_END 1330 1331 #undef __FUNC__ 1332 #define __FUNC__ "MatSeqSBAIJSetColumnIndices" 1333 /*@ 1334 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1335 in the matrix. 1336 1337 Input Parameters: 1338 + mat - the SeqSBAIJ matrix 1339 - indices - the column indices 1340 1341 Level: advanced 1342 1343 Notes: 1344 This can be called if you have precomputed the nonzero structure of the 1345 matrix and want to provide it to the matrix object to improve the performance 1346 of the MatSetValues() operation. 1347 1348 You MUST have set the correct numbers of nonzeros per row in the call to 1349 MatCreateSeqSBAIJ(). 1350 1351 MUST be called before any calls to MatSetValues() 1352 1353 .seealso: MatCreateSeqSBAIJ 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,*s_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 s_browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(s_browlengths); 1840 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1841 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1842 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1843 masked = mask + mbs; 1844 rowcount = 0; nzcount = 0; 1845 for (i=0; i<mbs; i++) { 1846 nmask = 0; 1847 for (j=0; j<bs; j++) { 1848 kmax = rowlengths[rowcount]; 1849 for (k=0; k<kmax; k++) { 1850 tmp = jj[nzcount++]/bs; /* block col. index */ 1851 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1852 } 1853 rowcount++; 1854 } 1855 s_browlengths[i] += nmask; 1856 browlengths[i] = 2*s_browlengths[i]; 1857 1858 /* zero out the mask elements we set */ 1859 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1860 } 1861 1862 /* create our matrix */ 1863 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1864 B = *A; 1865 a = (Mat_SeqSBAIJ*)B->data; 1866 1867 /* set matrix "i" values */ 1868 a->i[0] = 0; 1869 for (i=1; i<= mbs; i++) { 1870 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1871 a->ilen[i-1] = s_browlengths[i-1]; 1872 } 1873 a->s_nz = 0; 1874 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1875 1876 /* read in nonzero values */ 1877 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1878 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1879 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1880 1881 /* set "a" and "j" values into matrix */ 1882 nzcount = 0; jcount = 0; 1883 for (i=0; i<mbs; i++) { 1884 nzcountb = nzcount; 1885 nmask = 0; 1886 for (j=0; j<bs; j++) { 1887 kmax = rowlengths[i*bs+j]; 1888 for (k=0; k<kmax; k++) { 1889 tmp = jj[nzcount++]/bs; /* block col. index */ 1890 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1891 } 1892 rowcount++; 1893 } 1894 /* sort the masked values */ 1895 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1896 1897 /* set "j" values into matrix */ 1898 maskcount = 1; 1899 for (j=0; j<nmask; j++) { 1900 a->j[jcount++] = masked[j]; 1901 mask[masked[j]] = maskcount++; 1902 } 1903 1904 /* set "a" values into matrix */ 1905 ishift = bs2*a->i[i]; 1906 for (j=0; j<bs; j++) { 1907 kmax = rowlengths[i*bs+j]; 1908 for (k=0; k<kmax; k++) { 1909 tmp = jj[nzcountb]/bs ; /* block col. index */ 1910 if (tmp >= i){ 1911 block = mask[tmp] - 1; 1912 point = jj[nzcountb] - bs*tmp; 1913 idx = ishift + bs2*block + j + bs*point; 1914 a->a[idx] = aa[nzcountb]; 1915 } 1916 nzcountb++; 1917 } 1918 } 1919 /* zero out the mask elements we set */ 1920 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1921 } 1922 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1923 1924 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1925 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1926 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1927 ierr = PetscFree(aa);CHKERRQ(ierr); 1928 ierr = PetscFree(jj);CHKERRQ(ierr); 1929 ierr = PetscFree(mask);CHKERRQ(ierr); 1930 1931 B->assembled = PETSC_TRUE; 1932 ierr = MatView_Private(B);CHKERRQ(ierr); 1933 PetscFunctionReturn(0); 1934 } 1935 1936 1937 1938 1939