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