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