1 /*$Id: baij.c,v 1.210 2000/07/10 03:39:42 bsmith Exp bsmith $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "petscsys.h" 8 #include "src/mat/impls/baij/seq/baij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 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 = ISGetLocalSize(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,MatLUInfo*,Mat*); 1145 EXTERN int MatLUFactor_SeqBAIJ(Mat,IS,IS,MatLUInfo*); 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->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1231 inA->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering; 1232 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering; 1233 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=1\n"); 1234 case 2: 1235 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 1236 inA->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 1237 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 1238 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1239 break; 1240 case 3: 1241 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 1242 inA->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; 1243 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; 1244 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1245 break; 1246 case 4: 1247 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1248 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1249 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1250 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1251 break; 1252 case 5: 1253 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1254 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1255 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1256 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1257 break; 1258 case 6: 1259 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1260 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1261 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1262 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1263 break; 1264 case 7: 1265 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1266 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1267 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1268 PLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1269 break; 1270 default: 1271 a->row = row; 1272 a->col = col; 1273 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1274 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1275 1276 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1277 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1278 PLogObjectParent(inA,a->icol); 1279 1280 if (!a->solve_work) { 1281 a->solve_work = (Scalar*)PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1282 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1283 } 1284 } 1285 1286 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1287 1288 PetscFunctionReturn(0); 1289 } 1290 #undef __FUNC__ 1291 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp_SeqBAIJ" 1292 int MatPrintHelp_SeqBAIJ(Mat A) 1293 { 1294 static PetscTruth called = PETSC_FALSE; 1295 MPI_Comm comm = A->comm; 1296 int ierr; 1297 1298 PetscFunctionBegin; 1299 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1300 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1301 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 EXTERN_C_BEGIN 1306 #undef __FUNC__ 1307 #define __FUNC__ /*<a name=""></a>*/"MatSeqBAIJSetColumnIndices_SeqBAIJ" 1308 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1309 { 1310 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1311 int i,nz,n; 1312 1313 PetscFunctionBegin; 1314 nz = baij->maxnz; 1315 n = baij->n; 1316 for (i=0; i<nz; i++) { 1317 baij->j[i] = indices[i]; 1318 } 1319 baij->nz = nz; 1320 for (i=0; i<n; i++) { 1321 baij->ilen[i] = baij->imax[i]; 1322 } 1323 1324 PetscFunctionReturn(0); 1325 } 1326 EXTERN_C_END 1327 1328 #undef __FUNC__ 1329 #define __FUNC__ /*<a name=""></a>*/"MatSeqBAIJSetColumnIndices" 1330 /*@ 1331 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1332 in the matrix. 1333 1334 Input Parameters: 1335 + mat - the SeqBAIJ matrix 1336 - indices - the column indices 1337 1338 Level: advanced 1339 1340 Notes: 1341 This can be called if you have precomputed the nonzero structure of the 1342 matrix and want to provide it to the matrix object to improve the performance 1343 of the MatSetValues() operation. 1344 1345 You MUST have set the correct numbers of nonzeros per row in the call to 1346 MatCreateSeqBAIJ(). 1347 1348 MUST be called before any calls to MatSetValues(); 1349 1350 @*/ 1351 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1352 { 1353 int ierr,(*f)(Mat,int *); 1354 1355 PetscFunctionBegin; 1356 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1357 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1358 if (f) { 1359 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1360 } else { 1361 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1362 } 1363 PetscFunctionReturn(0); 1364 } 1365 1366 /* -------------------------------------------------------------------*/ 1367 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1368 MatGetRow_SeqBAIJ, 1369 MatRestoreRow_SeqBAIJ, 1370 MatMult_SeqBAIJ_N, 1371 MatMultAdd_SeqBAIJ_N, 1372 MatMultTranspose_SeqBAIJ, 1373 MatMultTransposeAdd_SeqBAIJ, 1374 MatSolve_SeqBAIJ_N, 1375 0, 1376 0, 1377 0, 1378 MatLUFactor_SeqBAIJ, 1379 0, 1380 0, 1381 MatTranspose_SeqBAIJ, 1382 MatGetInfo_SeqBAIJ, 1383 MatEqual_SeqBAIJ, 1384 MatGetDiagonal_SeqBAIJ, 1385 MatDiagonalScale_SeqBAIJ, 1386 MatNorm_SeqBAIJ, 1387 0, 1388 MatAssemblyEnd_SeqBAIJ, 1389 0, 1390 MatSetOption_SeqBAIJ, 1391 MatZeroEntries_SeqBAIJ, 1392 MatZeroRows_SeqBAIJ, 1393 MatLUFactorSymbolic_SeqBAIJ, 1394 MatLUFactorNumeric_SeqBAIJ_N, 1395 0, 1396 0, 1397 MatGetSize_SeqBAIJ, 1398 MatGetSize_SeqBAIJ, 1399 MatGetOwnershipRange_SeqBAIJ, 1400 MatILUFactorSymbolic_SeqBAIJ, 1401 0, 1402 0, 1403 0, 1404 MatDuplicate_SeqBAIJ, 1405 0, 1406 0, 1407 MatILUFactor_SeqBAIJ, 1408 0, 1409 0, 1410 MatGetSubMatrices_SeqBAIJ, 1411 MatIncreaseOverlap_SeqBAIJ, 1412 MatGetValues_SeqBAIJ, 1413 0, 1414 MatPrintHelp_SeqBAIJ, 1415 MatScale_SeqBAIJ, 1416 0, 1417 0, 1418 0, 1419 MatGetBlockSize_SeqBAIJ, 1420 MatGetRowIJ_SeqBAIJ, 1421 MatRestoreRowIJ_SeqBAIJ, 1422 0, 1423 0, 1424 0, 1425 0, 1426 0, 1427 0, 1428 MatSetValuesBlocked_SeqBAIJ, 1429 MatGetSubMatrix_SeqBAIJ, 1430 MatDestroy_SeqBAIJ, 1431 MatView_SeqBAIJ, 1432 MatGetMaps_Petsc}; 1433 1434 EXTERN_C_BEGIN 1435 #undef __FUNC__ 1436 #define __FUNC__ /*<a name=""></a>*/"MatStoreValues_SeqBAIJ" 1437 int MatStoreValues_SeqBAIJ(Mat mat) 1438 { 1439 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1440 int nz = aij->i[aij->m]*aij->bs*aij->bs2; 1441 int ierr; 1442 1443 PetscFunctionBegin; 1444 if (aij->nonew != 1) { 1445 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1446 } 1447 1448 /* allocate space for values if not already there */ 1449 if (!aij->saved_values) { 1450 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1451 } 1452 1453 /* copy values over */ 1454 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1455 PetscFunctionReturn(0); 1456 } 1457 EXTERN_C_END 1458 1459 EXTERN_C_BEGIN 1460 #undef __FUNC__ 1461 #define __FUNC__ /*<a name=""></a>*/"MatRetrieveValues_SeqBAIJ" 1462 int MatRetrieveValues_SeqBAIJ(Mat mat) 1463 { 1464 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1465 int nz = aij->i[aij->m]*aij->bs*aij->bs2,ierr; 1466 1467 PetscFunctionBegin; 1468 if (aij->nonew != 1) { 1469 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1470 } 1471 if (!aij->saved_values) { 1472 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1473 } 1474 1475 /* copy values over */ 1476 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1477 PetscFunctionReturn(0); 1478 } 1479 EXTERN_C_END 1480 1481 #undef __FUNC__ 1482 #define __FUNC__ /*<a name=""></a>*/"MatCreateSeqBAIJ" 1483 /*@C 1484 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1485 compressed row) format. For good matrix assembly performance the 1486 user should preallocate the matrix storage by setting the parameter nz 1487 (or the array nnz). By setting these parameters accurately, performance 1488 during matrix assembly can be increased by more than a factor of 50. 1489 1490 Collective on MPI_Comm 1491 1492 Input Parameters: 1493 + comm - MPI communicator, set to PETSC_COMM_SELF 1494 . bs - size of block 1495 . m - number of rows 1496 . n - number of columns 1497 . nz - number of block nonzeros per block row (same for all rows) 1498 - nnz - array containing the number of block nonzeros in the various block rows 1499 (possibly different for each block row) or PETSC_NULL 1500 1501 Output Parameter: 1502 . A - the matrix 1503 1504 Options Database Keys: 1505 . -mat_no_unroll - uses code that does not unroll the loops in the 1506 block calculations (much slower) 1507 . -mat_block_size - size of the blocks to use 1508 1509 Level: intermediate 1510 1511 Notes: 1512 The block AIJ format is fully compatible with standard Fortran 77 1513 storage. That is, the stored row and column indices can begin at 1514 either one (as in Fortran) or zero. See the users' manual for details. 1515 1516 Specify the preallocated storage with either nz or nnz (not both). 1517 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1518 allocation. For additional details, see the users manual chapter on 1519 matrices. 1520 1521 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1522 @*/ 1523 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1524 { 1525 Mat B; 1526 Mat_SeqBAIJ *b; 1527 int i,len,ierr,mbs,nbs,bs2,size,newbs = bs; 1528 PetscTruth flg; 1529 1530 PetscFunctionBegin; 1531 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1532 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1533 1534 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&newbs,PETSC_NULL);CHKERRQ(ierr); 1535 if (nnz && newbs != bs) { 1536 SETERRQ(1,1,"Cannot change blocksize from command line if setting nnz"); 1537 } 1538 1539 mbs = m/bs; 1540 nbs = n/bs; 1541 bs2 = bs*bs; 1542 1543 if (mbs*bs!=m || nbs*bs!=n) { 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,MATSEQBAIJ,"Mat",comm,MatDestroy,MatView); 1556 PLogObjectCreate(B); 1557 B->data = (void*)(b = PetscNew(Mat_SeqBAIJ));CHKPTRQ(b); 1558 ierr = PetscMemzero(b,sizeof(Mat_SeqBAIJ));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->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1565 B->ops->solve = MatSolve_SeqBAIJ_1; 1566 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 1567 B->ops->mult = MatMult_SeqBAIJ_1; 1568 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1569 break; 1570 case 2: 1571 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1572 B->ops->solve = MatSolve_SeqBAIJ_2; 1573 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 1574 B->ops->mult = MatMult_SeqBAIJ_2; 1575 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1576 break; 1577 case 3: 1578 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1579 B->ops->solve = MatSolve_SeqBAIJ_3; 1580 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3; 1581 B->ops->mult = MatMult_SeqBAIJ_3; 1582 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1583 break; 1584 case 4: 1585 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1586 B->ops->solve = MatSolve_SeqBAIJ_4; 1587 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4; 1588 B->ops->mult = MatMult_SeqBAIJ_4; 1589 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1590 break; 1591 case 5: 1592 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1593 B->ops->solve = MatSolve_SeqBAIJ_5; 1594 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5; 1595 B->ops->mult = MatMult_SeqBAIJ_5; 1596 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1597 break; 1598 case 6: 1599 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 1600 B->ops->solve = MatSolve_SeqBAIJ_6; 1601 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6; 1602 B->ops->mult = MatMult_SeqBAIJ_6; 1603 B->ops->multadd = MatMultAdd_SeqBAIJ_6; 1604 break; 1605 case 7: 1606 B->ops->mult = MatMult_SeqBAIJ_7; 1607 B->ops->solve = MatSolve_SeqBAIJ_7; 1608 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7; 1609 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1610 break; 1611 } 1612 } 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->factor = A->factor; 1712 c->row = 0; 1713 c->col = 0; 1714 c->icol = 0; 1715 c->saved_values = 0; 1716 c->keepzeroedrows = a->keepzeroedrows; 1717 C->assembled = PETSC_TRUE; 1718 1719 c->m = C->m = a->m; 1720 c->n = C->n = a->n; 1721 C->M = a->m; 1722 C->N = a->n; 1723 1724 c->bs = a->bs; 1725 c->bs2 = a->bs2; 1726 c->mbs = a->mbs; 1727 c->nbs = a->nbs; 1728 1729 c->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->imax); 1730 c->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->ilen); 1731 for (i=0; i<mbs; i++) { 1732 c->imax[i] = a->imax[i]; 1733 c->ilen[i] = a->ilen[i]; 1734 } 1735 1736 /* allocate the matrix space */ 1737 c->singlemalloc = PETSC_TRUE; 1738 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1739 c->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(c->a); 1740 c->j = (int*)(c->a + nz*bs2); 1741 c->i = c->j + nz; 1742 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1743 if (mbs > 0) { 1744 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1745 if (cpvalues == MAT_COPY_VALUES) { 1746 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1747 } else { 1748 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1749 } 1750 } 1751 1752 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1753 c->sorted = a->sorted; 1754 c->roworiented = a->roworiented; 1755 c->nonew = a->nonew; 1756 1757 if (a->diag) { 1758 c->diag = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->diag); 1759 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1760 for (i=0; i<mbs; i++) { 1761 c->diag[i] = a->diag[i]; 1762 } 1763 } else c->diag = 0; 1764 c->nz = a->nz; 1765 c->maxnz = a->maxnz; 1766 c->solve_work = 0; 1767 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1768 c->mult_work = 0; 1769 *B = C; 1770 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1771 PetscFunctionReturn(0); 1772 } 1773 1774 #undef __FUNC__ 1775 #define __FUNC__ /*<a name=""></a>*/"MatLoad_SeqBAIJ" 1776 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1777 { 1778 Mat_SeqBAIJ *a; 1779 Mat B; 1780 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1781 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1782 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1783 int *masked,nmask,tmp,bs2,ishift; 1784 Scalar *aa; 1785 MPI_Comm comm = ((PetscObject)viewer)->comm; 1786 1787 PetscFunctionBegin; 1788 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1789 bs2 = bs*bs; 1790 1791 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1792 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1793 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1794 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1795 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1796 M = header[1]; N = header[2]; nz = header[3]; 1797 1798 if (header[3] < 0) { 1799 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1800 } 1801 1802 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1803 1804 /* 1805 This code adds extra rows to make sure the number of rows is 1806 divisible by the blocksize 1807 */ 1808 mbs = M/bs; 1809 extra_rows = bs - M + bs*(mbs); 1810 if (extra_rows == bs) extra_rows = 0; 1811 else mbs++; 1812 if (extra_rows) { 1813 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1814 } 1815 1816 /* read in row lengths */ 1817 rowlengths = (int*)PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1818 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1819 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1820 1821 /* read in column indices */ 1822 jj = (int*)PetscMalloc((nz+extra_rows)*sizeof(int));CHKPTRQ(jj); 1823 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1824 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1825 1826 /* loop over row lengths determining block row lengths */ 1827 browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1828 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1829 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1830 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1831 masked = mask + mbs; 1832 rowcount = 0; nzcount = 0; 1833 for (i=0; i<mbs; i++) { 1834 nmask = 0; 1835 for (j=0; j<bs; j++) { 1836 kmax = rowlengths[rowcount]; 1837 for (k=0; k<kmax; k++) { 1838 tmp = jj[nzcount++]/bs; 1839 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1840 } 1841 rowcount++; 1842 } 1843 browlengths[i] += nmask; 1844 /* zero out the mask elements we set */ 1845 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1846 } 1847 1848 /* create our matrix */ 1849 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1850 B = *A; 1851 a = (Mat_SeqBAIJ*)B->data; 1852 1853 /* set matrix "i" values */ 1854 a->i[0] = 0; 1855 for (i=1; i<= mbs; i++) { 1856 a->i[i] = a->i[i-1] + browlengths[i-1]; 1857 a->ilen[i-1] = browlengths[i-1]; 1858 } 1859 a->nz = 0; 1860 for (i=0; i<mbs; i++) a->nz += browlengths[i]; 1861 1862 /* read in nonzero values */ 1863 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1864 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1865 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1866 1867 /* set "a" and "j" values into matrix */ 1868 nzcount = 0; jcount = 0; 1869 for (i=0; i<mbs; i++) { 1870 nzcountb = nzcount; 1871 nmask = 0; 1872 for (j=0; j<bs; j++) { 1873 kmax = rowlengths[i*bs+j]; 1874 for (k=0; k<kmax; k++) { 1875 tmp = jj[nzcount++]/bs; 1876 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1877 } 1878 rowcount++; 1879 } 1880 /* sort the masked values */ 1881 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1882 1883 /* set "j" values into matrix */ 1884 maskcount = 1; 1885 for (j=0; j<nmask; j++) { 1886 a->j[jcount++] = masked[j]; 1887 mask[masked[j]] = maskcount++; 1888 } 1889 /* set "a" values into matrix */ 1890 ishift = bs2*a->i[i]; 1891 for (j=0; j<bs; j++) { 1892 kmax = rowlengths[i*bs+j]; 1893 for (k=0; k<kmax; k++) { 1894 tmp = jj[nzcountb]/bs ; 1895 block = mask[tmp] - 1; 1896 point = jj[nzcountb] - bs*tmp; 1897 idx = ishift + bs2*block + j + bs*point; 1898 a->a[idx] = (MatScalar)aa[nzcountb++]; 1899 } 1900 } 1901 /* zero out the mask elements we set */ 1902 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1903 } 1904 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1905 1906 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1907 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1908 ierr = PetscFree(aa);CHKERRQ(ierr); 1909 ierr = PetscFree(jj);CHKERRQ(ierr); 1910 ierr = PetscFree(mask);CHKERRQ(ierr); 1911 1912 B->assembled = PETSC_TRUE; 1913 1914 ierr = MatView_Private(B);CHKERRQ(ierr); 1915 PetscFunctionReturn(0); 1916 } 1917 1918 1919 1920