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