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