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