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