1 /*$Id: baij.c,v 1.228 2001/06/21 22:55:30 buschelm 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 #if defined(PETSC_HAVE_SSE) 1203 PetscTruth sse_enabled; 1204 ierr = PetscSSEIsEnabled(&sse_enabled);CHKERRQ(ierr); 1205 if (sse_enabled) { 1206 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_SSE; 1207 } else { 1208 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1209 } 1210 #else 1211 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1212 #endif 1213 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1214 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4_NaturalOrdering; 1215 PetscLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1216 break; 1217 case 5: 1218 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1219 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1220 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; 1221 PetscLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1222 break; 1223 case 6: 1224 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 1225 inA->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 1226 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 1227 PetscLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1228 break; 1229 case 7: 1230 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 1231 inA->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1232 inA->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1233 PetscLogInfo(inA,"MatILUFactor_SeqBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1234 break; 1235 default: 1236 a->row = row; 1237 a->col = col; 1238 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1239 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1240 1241 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1242 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1243 PetscLogObjectParent(inA,a->icol); 1244 1245 if (!a->solve_work) { 1246 ierr = PetscMalloc((inA->m+a->bs)*sizeof(Scalar),&a->solve_work);CHKERRQ(ierr); 1247 PetscLogObjectMemory(inA,(inA->m+a->bs)*sizeof(Scalar)); 1248 } 1249 } 1250 1251 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1252 1253 PetscFunctionReturn(0); 1254 } 1255 #undef __FUNCT__ 1256 #define __FUNCT__ "MatPrintHelp_SeqBAIJ" 1257 int MatPrintHelp_SeqBAIJ(Mat A) 1258 { 1259 static PetscTruth called = PETSC_FALSE; 1260 MPI_Comm comm = A->comm; 1261 int ierr; 1262 1263 PetscFunctionBegin; 1264 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1265 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1266 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 1270 EXTERN_C_BEGIN 1271 #undef __FUNCT__ 1272 #define __FUNCT__ "MatSeqBAIJSetColumnIndices_SeqBAIJ" 1273 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1274 { 1275 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1276 int i,nz,n; 1277 1278 PetscFunctionBegin; 1279 nz = baij->maxnz; 1280 n = mat->n; 1281 for (i=0; i<nz; i++) { 1282 baij->j[i] = indices[i]; 1283 } 1284 baij->nz = nz; 1285 for (i=0; i<n; i++) { 1286 baij->ilen[i] = baij->imax[i]; 1287 } 1288 1289 PetscFunctionReturn(0); 1290 } 1291 EXTERN_C_END 1292 1293 #undef __FUNCT__ 1294 #define __FUNCT__ "MatSeqBAIJSetColumnIndices" 1295 /*@ 1296 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1297 in the matrix. 1298 1299 Input Parameters: 1300 + mat - the SeqBAIJ matrix 1301 - indices - the column indices 1302 1303 Level: advanced 1304 1305 Notes: 1306 This can be called if you have precomputed the nonzero structure of the 1307 matrix and want to provide it to the matrix object to improve the performance 1308 of the MatSetValues() operation. 1309 1310 You MUST have set the correct numbers of nonzeros per row in the call to 1311 MatCreateSeqBAIJ(). 1312 1313 MUST be called before any calls to MatSetValues(); 1314 1315 @*/ 1316 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1317 { 1318 int ierr,(*f)(Mat,int *); 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1322 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void (**)())&f);CHKERRQ(ierr); 1323 if (f) { 1324 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1325 } else { 1326 SETERRQ(1,"Wrong type of matrix to set column indices"); 1327 } 1328 PetscFunctionReturn(0); 1329 } 1330 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "MatGetRowMax_SeqBAIJ" 1333 int MatGetRowMax_SeqBAIJ(Mat A,Vec v) 1334 { 1335 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1336 int ierr,i,j,n,row,bs,*ai,*aj,mbs; 1337 PetscReal atmp; 1338 Scalar *x,zero = 0.0; 1339 MatScalar *aa; 1340 int ncols,brow,krow,kcol; 1341 1342 PetscFunctionBegin; 1343 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1344 bs = a->bs; 1345 aa = a->a; 1346 ai = a->i; 1347 aj = a->j; 1348 mbs = a->mbs; 1349 1350 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1351 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1352 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1353 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1354 for (i=0; i<mbs; i++) { 1355 ncols = ai[1] - ai[0]; ai++; 1356 brow = bs*i; 1357 for (j=0; j<ncols; j++){ 1358 /* bcol = bs*(*aj); */ 1359 for (kcol=0; kcol<bs; kcol++){ 1360 for (krow=0; krow<bs; krow++){ 1361 atmp = PetscAbsScalar(*aa); aa++; 1362 row = brow + krow; /* row index */ 1363 /* printf("val[%d,%d]: %g\n",row,bcol+kcol,atmp); */ 1364 if (PetscAbsScalar(x[row]) < atmp) x[row] = atmp; 1365 } 1366 } 1367 aj++; 1368 } 1369 } 1370 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1371 PetscFunctionReturn(0); 1372 } 1373 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "MatSetUpPreallocation_SeqBAIJ" 1376 int MatSetUpPreallocation_SeqBAIJ(Mat A) 1377 { 1378 int ierr; 1379 1380 PetscFunctionBegin; 1381 ierr = MatSeqBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 #undef __FUNCT__ 1386 #define __FUNCT__ "MatGetArray_SeqBAIJ" 1387 int MatGetArray_SeqBAIJ(Mat A,Scalar **array) 1388 { 1389 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1390 PetscFunctionBegin; 1391 *array = a->a; 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "MatRestoreArray_SeqBAIJ" 1397 int MatRestoreArray_SeqBAIJ(Mat A,Scalar **array) 1398 { 1399 PetscFunctionBegin; 1400 PetscFunctionReturn(0); 1401 } 1402 1403 /* -------------------------------------------------------------------*/ 1404 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1405 MatGetRow_SeqBAIJ, 1406 MatRestoreRow_SeqBAIJ, 1407 MatMult_SeqBAIJ_N, 1408 MatMultAdd_SeqBAIJ_N, 1409 MatMultTranspose_SeqBAIJ, 1410 MatMultTransposeAdd_SeqBAIJ, 1411 MatSolve_SeqBAIJ_N, 1412 0, 1413 0, 1414 0, 1415 MatLUFactor_SeqBAIJ, 1416 0, 1417 0, 1418 MatTranspose_SeqBAIJ, 1419 MatGetInfo_SeqBAIJ, 1420 MatEqual_SeqBAIJ, 1421 MatGetDiagonal_SeqBAIJ, 1422 MatDiagonalScale_SeqBAIJ, 1423 MatNorm_SeqBAIJ, 1424 0, 1425 MatAssemblyEnd_SeqBAIJ, 1426 0, 1427 MatSetOption_SeqBAIJ, 1428 MatZeroEntries_SeqBAIJ, 1429 MatZeroRows_SeqBAIJ, 1430 MatLUFactorSymbolic_SeqBAIJ, 1431 MatLUFactorNumeric_SeqBAIJ_N, 1432 0, 1433 0, 1434 MatSetUpPreallocation_SeqBAIJ, 1435 0, 1436 MatGetOwnershipRange_SeqBAIJ, 1437 MatILUFactorSymbolic_SeqBAIJ, 1438 0, 1439 MatGetArray_SeqBAIJ, 1440 MatRestoreArray_SeqBAIJ, 1441 MatDuplicate_SeqBAIJ, 1442 0, 1443 0, 1444 MatILUFactor_SeqBAIJ, 1445 0, 1446 0, 1447 MatGetSubMatrices_SeqBAIJ, 1448 MatIncreaseOverlap_SeqBAIJ, 1449 MatGetValues_SeqBAIJ, 1450 0, 1451 MatPrintHelp_SeqBAIJ, 1452 MatScale_SeqBAIJ, 1453 0, 1454 0, 1455 0, 1456 MatGetBlockSize_SeqBAIJ, 1457 MatGetRowIJ_SeqBAIJ, 1458 MatRestoreRowIJ_SeqBAIJ, 1459 0, 1460 0, 1461 0, 1462 0, 1463 0, 1464 0, 1465 MatSetValuesBlocked_SeqBAIJ, 1466 MatGetSubMatrix_SeqBAIJ, 1467 MatDestroy_SeqBAIJ, 1468 MatView_SeqBAIJ, 1469 MatGetMaps_Petsc, 1470 0, 1471 0, 1472 0, 1473 0, 1474 0, 1475 0, 1476 MatGetRowMax_SeqBAIJ, 1477 MatConvert_Basic}; 1478 1479 EXTERN_C_BEGIN 1480 #undef __FUNCT__ 1481 #define __FUNCT__ "MatStoreValues_SeqBAIJ" 1482 int MatStoreValues_SeqBAIJ(Mat mat) 1483 { 1484 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1485 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1486 int ierr; 1487 1488 PetscFunctionBegin; 1489 if (aij->nonew != 1) { 1490 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1491 } 1492 1493 /* allocate space for values if not already there */ 1494 if (!aij->saved_values) { 1495 ierr = PetscMalloc((nz+1)*sizeof(Scalar),&aij->saved_values);CHKERRQ(ierr); 1496 } 1497 1498 /* copy values over */ 1499 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1500 PetscFunctionReturn(0); 1501 } 1502 EXTERN_C_END 1503 1504 EXTERN_C_BEGIN 1505 #undef __FUNCT__ 1506 #define __FUNCT__ "MatRetrieveValues_SeqBAIJ" 1507 int MatRetrieveValues_SeqBAIJ(Mat mat) 1508 { 1509 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1510 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1511 1512 PetscFunctionBegin; 1513 if (aij->nonew != 1) { 1514 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1515 } 1516 if (!aij->saved_values) { 1517 SETERRQ(1,"Must call MatStoreValues(A);first"); 1518 } 1519 1520 /* copy values over */ 1521 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1522 PetscFunctionReturn(0); 1523 } 1524 EXTERN_C_END 1525 1526 EXTERN_C_BEGIN 1527 extern int MatConvert_SeqBAIJ_SeqAIJ(Mat,MatType,Mat*); 1528 EXTERN_C_END 1529 1530 EXTERN_C_BEGIN 1531 #undef __FUNCT__ 1532 #define __FUNCT__ "MatCreate_SeqBAIJ" 1533 int MatCreate_SeqBAIJ(Mat B) 1534 { 1535 int ierr,size; 1536 Mat_SeqBAIJ *b; 1537 1538 PetscFunctionBegin; 1539 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1540 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1541 1542 B->m = B->M = PetscMax(B->m,B->M); 1543 B->n = B->N = PetscMax(B->n,B->N); 1544 ierr = PetscNew(Mat_SeqBAIJ,&b);CHKERRQ(ierr); 1545 B->data = (void*)b; 1546 ierr = PetscMemzero(b,sizeof(Mat_SeqBAIJ));CHKERRQ(ierr); 1547 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1548 B->factor = 0; 1549 B->lupivotthreshold = 1.0; 1550 B->mapping = 0; 1551 b->row = 0; 1552 b->col = 0; 1553 b->icol = 0; 1554 b->reallocs = 0; 1555 b->saved_values = 0; 1556 #if defined(PETSC_USE_MAT_SINGLE) 1557 b->setvalueslen = 0; 1558 b->setvaluescopy = PETSC_NULL; 1559 #endif 1560 1561 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1562 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1563 1564 b->sorted = PETSC_FALSE; 1565 b->roworiented = PETSC_TRUE; 1566 b->nonew = 0; 1567 b->diag = 0; 1568 b->solve_work = 0; 1569 b->mult_work = 0; 1570 b->spptr = 0; 1571 B->info.nz_unneeded = (PetscReal)b->maxnz; 1572 b->keepzeroedrows = PETSC_FALSE; 1573 1574 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1575 "MatStoreValues_SeqBAIJ", 1576 MatStoreValues_SeqBAIJ);CHKERRQ(ierr); 1577 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1578 "MatRetrieveValues_SeqBAIJ", 1579 MatRetrieveValues_SeqBAIJ);CHKERRQ(ierr); 1580 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1581 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1582 MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1583 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_SeqBAIJ_SeqAIJ_C", 1584 "MatConvert_SeqBAIJ_SeqAIJ", 1585 MatConvert_SeqBAIJ_SeqAIJ);CHKERRQ(ierr); 1586 PetscFunctionReturn(0); 1587 } 1588 EXTERN_C_END 1589 1590 #undef __FUNCT__ 1591 #define __FUNCT__ "MatDuplicate_SeqBAIJ" 1592 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1593 { 1594 Mat C; 1595 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ*)A->data; 1596 int i,len,mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1597 1598 PetscFunctionBegin; 1599 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1600 1601 *B = 0; 1602 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1603 ierr = MatSetType(C,MATSEQBAIJ);CHKERRQ(ierr); 1604 c = (Mat_SeqBAIJ*)C->data; 1605 1606 c->bs = a->bs; 1607 c->bs2 = a->bs2; 1608 c->mbs = a->mbs; 1609 c->nbs = a->nbs; 1610 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1611 1612 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 1613 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 1614 for (i=0; i<mbs; i++) { 1615 c->imax[i] = a->imax[i]; 1616 c->ilen[i] = a->ilen[i]; 1617 } 1618 1619 /* allocate the matrix space */ 1620 c->singlemalloc = PETSC_TRUE; 1621 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1622 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 1623 c->j = (int*)(c->a + nz*bs2); 1624 c->i = c->j + nz; 1625 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1626 if (mbs > 0) { 1627 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1628 if (cpvalues == MAT_COPY_VALUES) { 1629 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1630 } else { 1631 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1632 } 1633 } 1634 1635 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1636 c->sorted = a->sorted; 1637 c->roworiented = a->roworiented; 1638 c->nonew = a->nonew; 1639 1640 if (a->diag) { 1641 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 1642 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1643 for (i=0; i<mbs; i++) { 1644 c->diag[i] = a->diag[i]; 1645 } 1646 } else c->diag = 0; 1647 c->nz = a->nz; 1648 c->maxnz = a->maxnz; 1649 c->solve_work = 0; 1650 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1651 c->mult_work = 0; 1652 C->preallocated = PETSC_TRUE; 1653 C->assembled = PETSC_TRUE; 1654 *B = C; 1655 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 EXTERN_C_BEGIN 1660 #undef __FUNCT__ 1661 #define __FUNCT__ "MatLoad_SeqBAIJ" 1662 int MatLoad_SeqBAIJ(PetscViewer viewer,MatType type,Mat *A) 1663 { 1664 Mat_SeqBAIJ *a; 1665 Mat B; 1666 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1667 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1668 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1669 int *masked,nmask,tmp,bs2,ishift; 1670 Scalar *aa; 1671 MPI_Comm comm = ((PetscObject)viewer)->comm; 1672 1673 PetscFunctionBegin; 1674 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1675 bs2 = bs*bs; 1676 1677 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1678 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1679 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1680 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1681 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1682 M = header[1]; N = header[2]; nz = header[3]; 1683 1684 if (header[3] < 0) { 1685 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBAIJ"); 1686 } 1687 1688 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1689 1690 /* 1691 This code adds extra rows to make sure the number of rows is 1692 divisible by the blocksize 1693 */ 1694 mbs = M/bs; 1695 extra_rows = bs - M + bs*(mbs); 1696 if (extra_rows == bs) extra_rows = 0; 1697 else mbs++; 1698 if (extra_rows) { 1699 PetscLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1700 } 1701 1702 /* read in row lengths */ 1703 ierr = PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);CHKERRQ(ierr); 1704 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1705 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1706 1707 /* read in column indices */ 1708 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),&jj);CHKERRQ(ierr); 1709 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1710 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1711 1712 /* loop over row lengths determining block row lengths */ 1713 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 1714 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1715 ierr = PetscMalloc(2*mbs*sizeof(int),&mask);CHKERRQ(ierr); 1716 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1717 masked = mask + mbs; 1718 rowcount = 0; nzcount = 0; 1719 for (i=0; i<mbs; i++) { 1720 nmask = 0; 1721 for (j=0; j<bs; j++) { 1722 kmax = rowlengths[rowcount]; 1723 for (k=0; k<kmax; k++) { 1724 tmp = jj[nzcount++]/bs; 1725 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1726 } 1727 rowcount++; 1728 } 1729 browlengths[i] += nmask; 1730 /* zero out the mask elements we set */ 1731 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1732 } 1733 1734 /* create our matrix */ 1735 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1736 B = *A; 1737 a = (Mat_SeqBAIJ*)B->data; 1738 1739 /* set matrix "i" values */ 1740 a->i[0] = 0; 1741 for (i=1; i<= mbs; i++) { 1742 a->i[i] = a->i[i-1] + browlengths[i-1]; 1743 a->ilen[i-1] = browlengths[i-1]; 1744 } 1745 a->nz = 0; 1746 for (i=0; i<mbs; i++) a->nz += browlengths[i]; 1747 1748 /* read in nonzero values */ 1749 ierr = PetscMalloc((nz+extra_rows)*sizeof(Scalar),&aa);CHKERRQ(ierr); 1750 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1751 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1752 1753 /* set "a" and "j" values into matrix */ 1754 nzcount = 0; jcount = 0; 1755 for (i=0; i<mbs; i++) { 1756 nzcountb = nzcount; 1757 nmask = 0; 1758 for (j=0; j<bs; j++) { 1759 kmax = rowlengths[i*bs+j]; 1760 for (k=0; k<kmax; k++) { 1761 tmp = jj[nzcount++]/bs; 1762 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1763 } 1764 } 1765 /* sort the masked values */ 1766 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1767 1768 /* set "j" values into matrix */ 1769 maskcount = 1; 1770 for (j=0; j<nmask; j++) { 1771 a->j[jcount++] = masked[j]; 1772 mask[masked[j]] = maskcount++; 1773 } 1774 /* set "a" values into matrix */ 1775 ishift = bs2*a->i[i]; 1776 for (j=0; j<bs; j++) { 1777 kmax = rowlengths[i*bs+j]; 1778 for (k=0; k<kmax; k++) { 1779 tmp = jj[nzcountb]/bs ; 1780 block = mask[tmp] - 1; 1781 point = jj[nzcountb] - bs*tmp; 1782 idx = ishift + bs2*block + j + bs*point; 1783 a->a[idx] = (MatScalar)aa[nzcountb++]; 1784 } 1785 } 1786 /* zero out the mask elements we set */ 1787 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1788 } 1789 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1790 1791 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1792 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1793 ierr = PetscFree(aa);CHKERRQ(ierr); 1794 ierr = PetscFree(jj);CHKERRQ(ierr); 1795 ierr = PetscFree(mask);CHKERRQ(ierr); 1796 1797 B->assembled = PETSC_TRUE; 1798 1799 ierr = MatView_Private(B);CHKERRQ(ierr); 1800 PetscFunctionReturn(0); 1801 } 1802 EXTERN_C_END 1803 1804 #undef __FUNCT__ 1805 #define __FUNCT__ "MatCreateSeqBAIJ" 1806 /*@C 1807 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1808 compressed row) format. For good matrix assembly performance the 1809 user should preallocate the matrix storage by setting the parameter nz 1810 (or the array nnz). By setting these parameters accurately, performance 1811 during matrix assembly can be increased by more than a factor of 50. 1812 1813 Collective on MPI_Comm 1814 1815 Input Parameters: 1816 + comm - MPI communicator, set to PETSC_COMM_SELF 1817 . bs - size of block 1818 . m - number of rows 1819 . n - number of columns 1820 . nz - number of nonzero blocks per block row (same for all rows) 1821 - nnz - array containing the number of nonzero blocks in the various block rows 1822 (possibly different for each block row) or PETSC_NULL 1823 1824 Output Parameter: 1825 . A - the matrix 1826 1827 Options Database Keys: 1828 . -mat_no_unroll - uses code that does not unroll the loops in the 1829 block calculations (much slower) 1830 . -mat_block_size - size of the blocks to use 1831 1832 Level: intermediate 1833 1834 Notes: 1835 A nonzero block is any block that as 1 or more nonzeros in it 1836 1837 The block AIJ format is fully compatible with standard Fortran 77 1838 storage. That is, the stored row and column indices can begin at 1839 either one (as in Fortran) or zero. See the users' manual for details. 1840 1841 Specify the preallocated storage with either nz or nnz (not both). 1842 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1843 allocation. For additional details, see the users manual chapter on 1844 matrices. 1845 1846 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1847 @*/ 1848 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1849 { 1850 int ierr; 1851 1852 PetscFunctionBegin; 1853 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1854 ierr = MatSetType(*A,MATSEQBAIJ);CHKERRQ(ierr); 1855 ierr = MatSeqBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1856 PetscFunctionReturn(0); 1857 } 1858 1859 #undef __FUNCT__ 1860 #define __FUNCT__ "MatSeqBAIJSetPreallocation" 1861 /*@C 1862 MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros 1863 per row in the matrix. For good matrix assembly performance the 1864 user should preallocate the matrix storage by setting the parameter nz 1865 (or the array nnz). By setting these parameters accurately, performance 1866 during matrix assembly can be increased by more than a factor of 50. 1867 1868 Collective on MPI_Comm 1869 1870 Input Parameters: 1871 + A - the matrix 1872 . bs - size of block 1873 . nz - number of block nonzeros per block row (same for all rows) 1874 - nnz - array containing the number of block nonzeros in the various block rows 1875 (possibly different for each block row) or PETSC_NULL 1876 1877 Options Database Keys: 1878 . -mat_no_unroll - uses code that does not unroll the loops in the 1879 block calculations (much slower) 1880 . -mat_block_size - size of the blocks to use 1881 1882 Level: intermediate 1883 1884 Notes: 1885 The block AIJ format is fully compatible with standard Fortran 77 1886 storage. That is, the stored row and column indices can begin at 1887 either one (as in Fortran) or zero. See the users' manual for details. 1888 1889 Specify the preallocated storage with either nz or nnz (not both). 1890 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1891 allocation. For additional details, see the users manual chapter on 1892 matrices. 1893 1894 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1895 @*/ 1896 int MatSeqBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1897 { 1898 Mat_SeqBAIJ *b; 1899 int i,len,ierr,mbs,nbs,bs2,newbs = bs; 1900 PetscTruth flg; 1901 1902 PetscFunctionBegin; 1903 ierr = PetscTypeCompare((PetscObject)B,MATSEQBAIJ,&flg);CHKERRQ(ierr); 1904 if (!flg) PetscFunctionReturn(0); 1905 1906 B->preallocated = PETSC_TRUE; 1907 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&newbs,PETSC_NULL);CHKERRQ(ierr); 1908 if (nnz && newbs != bs) { 1909 SETERRQ(1,"Cannot change blocksize from command line if setting nnz"); 1910 } 1911 bs = newbs; 1912 1913 mbs = B->m/bs; 1914 nbs = B->n/bs; 1915 bs2 = bs*bs; 1916 1917 if (mbs*bs!=B->m || nbs*bs!=B->n) { 1918 SETERRQ3(PETSC_ERR_ARG_SIZ,"Number rows %d, cols %d must be divisible by blocksize %d",B->m,B->n,bs); 1919 } 1920 1921 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1922 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 1923 if (nnz) { 1924 for (i=0; i<mbs; i++) { 1925 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1926 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); 1927 } 1928 } 1929 1930 b = (Mat_SeqBAIJ*)B->data; 1931 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1932 if (!flg) { 1933 switch (bs) { 1934 case 1: 1935 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1936 B->ops->solve = MatSolve_SeqBAIJ_1; 1937 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 1938 B->ops->mult = MatMult_SeqBAIJ_1; 1939 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1940 break; 1941 case 2: 1942 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1943 B->ops->solve = MatSolve_SeqBAIJ_2; 1944 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 1945 B->ops->mult = MatMult_SeqBAIJ_2; 1946 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1947 break; 1948 case 3: 1949 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1950 B->ops->solve = MatSolve_SeqBAIJ_3; 1951 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3; 1952 B->ops->mult = MatMult_SeqBAIJ_3; 1953 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1954 break; 1955 case 4: 1956 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1957 B->ops->solve = MatSolve_SeqBAIJ_4; 1958 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4; 1959 B->ops->mult = MatMult_SeqBAIJ_4; 1960 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1961 break; 1962 case 5: 1963 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1964 B->ops->solve = MatSolve_SeqBAIJ_5; 1965 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5; 1966 B->ops->mult = MatMult_SeqBAIJ_5; 1967 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1968 break; 1969 case 6: 1970 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 1971 B->ops->solve = MatSolve_SeqBAIJ_6; 1972 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6; 1973 B->ops->mult = MatMult_SeqBAIJ_6; 1974 B->ops->multadd = MatMultAdd_SeqBAIJ_6; 1975 break; 1976 case 7: 1977 B->ops->mult = MatMult_SeqBAIJ_7; 1978 B->ops->solve = MatSolve_SeqBAIJ_7; 1979 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7; 1980 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1981 break; 1982 default: 1983 B->ops->mult = MatMult_SeqBAIJ_N; 1984 B->ops->solve = MatSolve_SeqBAIJ_N; 1985 B->ops->multadd = MatMultAdd_SeqBAIJ_N; 1986 break; 1987 } 1988 } 1989 b->bs = bs; 1990 b->mbs = mbs; 1991 b->nbs = nbs; 1992 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 1993 if (!nnz) { 1994 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1995 else if (nz <= 0) nz = 1; 1996 for (i=0; i<mbs; i++) b->imax[i] = nz; 1997 nz = nz*mbs; 1998 } else { 1999 nz = 0; 2000 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2001 } 2002 2003 /* allocate the matrix space */ 2004 len = nz*sizeof(int) + nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 2005 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2006 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2007 b->j = (int*)(b->a + nz*bs2); 2008 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2009 b->i = b->j + nz; 2010 b->singlemalloc = PETSC_TRUE; 2011 2012 b->i[0] = 0; 2013 for (i=1; i<mbs+1; i++) { 2014 b->i[i] = b->i[i-1] + b->imax[i-1]; 2015 } 2016 2017 /* b->ilen will count nonzeros in each block row so far. */ 2018 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2019 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 2020 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 2021 2022 b->bs = bs; 2023 b->bs2 = bs2; 2024 b->mbs = mbs; 2025 b->nz = 0; 2026 b->maxnz = nz*bs2; 2027 B->info.nz_unneeded = (PetscReal)b->maxnz; 2028 PetscFunctionReturn(0); 2029 } 2030 2031