1 /*$Id: baij.c,v 1.217 2001/01/19 23:20:39 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] + 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 /* -------------------------------------------------------------------*/ 1359 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1360 MatGetRow_SeqBAIJ, 1361 MatRestoreRow_SeqBAIJ, 1362 MatMult_SeqBAIJ_N, 1363 MatMultAdd_SeqBAIJ_N, 1364 MatMultTranspose_SeqBAIJ, 1365 MatMultTransposeAdd_SeqBAIJ, 1366 MatSolve_SeqBAIJ_N, 1367 0, 1368 0, 1369 0, 1370 MatLUFactor_SeqBAIJ, 1371 0, 1372 0, 1373 MatTranspose_SeqBAIJ, 1374 MatGetInfo_SeqBAIJ, 1375 MatEqual_SeqBAIJ, 1376 MatGetDiagonal_SeqBAIJ, 1377 MatDiagonalScale_SeqBAIJ, 1378 MatNorm_SeqBAIJ, 1379 0, 1380 MatAssemblyEnd_SeqBAIJ, 1381 0, 1382 MatSetOption_SeqBAIJ, 1383 MatZeroEntries_SeqBAIJ, 1384 MatZeroRows_SeqBAIJ, 1385 MatLUFactorSymbolic_SeqBAIJ, 1386 MatLUFactorNumeric_SeqBAIJ_N, 1387 0, 1388 0, 1389 MatSetUpPreallocation_SeqBAIJ, 1390 0, 1391 MatGetOwnershipRange_SeqBAIJ, 1392 MatILUFactorSymbolic_SeqBAIJ, 1393 0, 1394 0, 1395 0, 1396 MatDuplicate_SeqBAIJ, 1397 0, 1398 0, 1399 MatILUFactor_SeqBAIJ, 1400 0, 1401 0, 1402 MatGetSubMatrices_SeqBAIJ, 1403 MatIncreaseOverlap_SeqBAIJ, 1404 MatGetValues_SeqBAIJ, 1405 0, 1406 MatPrintHelp_SeqBAIJ, 1407 MatScale_SeqBAIJ, 1408 0, 1409 0, 1410 0, 1411 MatGetBlockSize_SeqBAIJ, 1412 MatGetRowIJ_SeqBAIJ, 1413 MatRestoreRowIJ_SeqBAIJ, 1414 0, 1415 0, 1416 0, 1417 0, 1418 0, 1419 0, 1420 MatSetValuesBlocked_SeqBAIJ, 1421 MatGetSubMatrix_SeqBAIJ, 1422 MatDestroy_SeqBAIJ, 1423 MatView_SeqBAIJ, 1424 MatGetMaps_Petsc, 1425 0, 1426 0, 1427 0, 1428 0, 1429 0, 1430 0, 1431 MatGetRowMax_SeqBAIJ, 1432 MatConvert_Basic}; 1433 1434 EXTERN_C_BEGIN 1435 #undef __FUNC__ 1436 #define __FUNC__ "MatStoreValues_SeqBAIJ" 1437 int MatStoreValues_SeqBAIJ(Mat mat) 1438 { 1439 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1440 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1441 int ierr; 1442 1443 PetscFunctionBegin; 1444 if (aij->nonew != 1) { 1445 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1446 } 1447 1448 /* allocate space for values if not already there */ 1449 if (!aij->saved_values) { 1450 ierr = PetscMalloc(nz*sizeof(Scalar),&aij->saved_values);CHKERRQ(ierr); 1451 } 1452 1453 /* copy values over */ 1454 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1455 PetscFunctionReturn(0); 1456 } 1457 EXTERN_C_END 1458 1459 EXTERN_C_BEGIN 1460 #undef __FUNC__ 1461 #define __FUNC__ "MatRetrieveValues_SeqBAIJ" 1462 int MatRetrieveValues_SeqBAIJ(Mat mat) 1463 { 1464 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1465 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1466 1467 PetscFunctionBegin; 1468 if (aij->nonew != 1) { 1469 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1470 } 1471 if (!aij->saved_values) { 1472 SETERRQ(1,"Must call MatStoreValues(A);first"); 1473 } 1474 1475 /* copy values over */ 1476 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1477 PetscFunctionReturn(0); 1478 } 1479 EXTERN_C_END 1480 1481 EXTERN_C_BEGIN 1482 extern int MatConvert_SeqBAIJ_SeqAIJ(Mat,MatType,Mat*); 1483 EXTERN_C_END 1484 1485 EXTERN_C_BEGIN 1486 #undef __FUNC__ 1487 #define __FUNC__ "MatCreate_SeqBAIJ" 1488 int MatCreate_SeqBAIJ(Mat B) 1489 { 1490 int ierr,size; 1491 Mat_SeqBAIJ *b; 1492 1493 PetscFunctionBegin; 1494 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1495 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1496 1497 B->m = B->M = PetscMax(B->m,B->M); 1498 B->n = B->N = PetscMax(B->n,B->N); 1499 ierr = PetscNew(Mat_SeqBAIJ,&b);CHKERRQ(ierr); 1500 B->data = (void*)b; 1501 ierr = PetscMemzero(b,sizeof(Mat_SeqBAIJ));CHKERRQ(ierr); 1502 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1503 B->factor = 0; 1504 B->lupivotthreshold = 1.0; 1505 B->mapping = 0; 1506 b->row = 0; 1507 b->col = 0; 1508 b->icol = 0; 1509 b->reallocs = 0; 1510 b->saved_values = 0; 1511 #if defined(PEYSC_USE_MAT_SINGLE) 1512 b->setvalueslen = 0; 1513 b->setvaluescopy = PETSC_NULL; 1514 #endif 1515 1516 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1517 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1518 1519 b->sorted = PETSC_FALSE; 1520 b->roworiented = PETSC_TRUE; 1521 b->nonew = 0; 1522 b->diag = 0; 1523 b->solve_work = 0; 1524 b->mult_work = 0; 1525 b->spptr = 0; 1526 B->info.nz_unneeded = (PetscReal)b->maxnz; 1527 b->keepzeroedrows = PETSC_FALSE; 1528 1529 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1530 "MatStoreValues_SeqBAIJ", 1531 MatStoreValues_SeqBAIJ);CHKERRQ(ierr); 1532 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1533 "MatRetrieveValues_SeqBAIJ", 1534 MatRetrieveValues_SeqBAIJ);CHKERRQ(ierr); 1535 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1536 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1537 MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1538 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_SeqBAIJ_SeqAIJ_C", 1539 "MatConvert_SeqBAIJ_SeqAIJ", 1540 MatConvert_SeqBAIJ_SeqAIJ);CHKERRQ(ierr); 1541 PetscFunctionReturn(0); 1542 } 1543 EXTERN_C_END 1544 1545 #undef __FUNC__ 1546 #define __FUNC__ "MatDuplicate_SeqBAIJ" 1547 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1548 { 1549 Mat C; 1550 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ*)A->data; 1551 int i,len,mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1552 1553 PetscFunctionBegin; 1554 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1555 1556 *B = 0; 1557 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1558 ierr = MatSetType(C,MATSEQBAIJ);CHKERRQ(ierr); 1559 c = (Mat_SeqBAIJ*)C->data; 1560 1561 c->bs = a->bs; 1562 c->bs2 = a->bs2; 1563 c->mbs = a->mbs; 1564 c->nbs = a->nbs; 1565 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1566 1567 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 1568 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 1569 for (i=0; i<mbs; i++) { 1570 c->imax[i] = a->imax[i]; 1571 c->ilen[i] = a->ilen[i]; 1572 } 1573 1574 /* allocate the matrix space */ 1575 c->singlemalloc = PETSC_TRUE; 1576 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1577 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 1578 c->j = (int*)(c->a + nz*bs2); 1579 c->i = c->j + nz; 1580 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1581 if (mbs > 0) { 1582 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1583 if (cpvalues == MAT_COPY_VALUES) { 1584 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1585 } else { 1586 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1587 } 1588 } 1589 1590 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1591 c->sorted = a->sorted; 1592 c->roworiented = a->roworiented; 1593 c->nonew = a->nonew; 1594 1595 if (a->diag) { 1596 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 1597 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1598 for (i=0; i<mbs; i++) { 1599 c->diag[i] = a->diag[i]; 1600 } 1601 } else c->diag = 0; 1602 c->nz = a->nz; 1603 c->maxnz = a->maxnz; 1604 c->solve_work = 0; 1605 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1606 c->mult_work = 0; 1607 C->preallocated = PETSC_TRUE; 1608 C->assembled = PETSC_TRUE; 1609 *B = C; 1610 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 EXTERN_C_BEGIN 1615 #undef __FUNC__ 1616 #define __FUNC__ "MatLoad_SeqBAIJ" 1617 int MatLoad_SeqBAIJ(PetscViewer viewer,MatType type,Mat *A) 1618 { 1619 Mat_SeqBAIJ *a; 1620 Mat B; 1621 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1622 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1623 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1624 int *masked,nmask,tmp,bs2,ishift; 1625 Scalar *aa; 1626 MPI_Comm comm = ((PetscObject)viewer)->comm; 1627 1628 PetscFunctionBegin; 1629 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1630 bs2 = bs*bs; 1631 1632 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1633 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1634 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1635 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1636 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1637 M = header[1]; N = header[2]; nz = header[3]; 1638 1639 if (header[3] < 0) { 1640 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBAIJ"); 1641 } 1642 1643 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1644 1645 /* 1646 This code adds extra rows to make sure the number of rows is 1647 divisible by the blocksize 1648 */ 1649 mbs = M/bs; 1650 extra_rows = bs - M + bs*(mbs); 1651 if (extra_rows == bs) extra_rows = 0; 1652 else mbs++; 1653 if (extra_rows) { 1654 PetscLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1655 } 1656 1657 /* read in row lengths */ 1658 ierr = PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);CHKERRQ(ierr); 1659 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1660 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1661 1662 /* read in column indices */ 1663 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),&jj);CHKERRQ(ierr); 1664 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1665 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1666 1667 /* loop over row lengths determining block row lengths */ 1668 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 1669 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1670 ierr = PetscMalloc(2*mbs*sizeof(int),&mask);CHKERRQ(ierr); 1671 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1672 masked = mask + mbs; 1673 rowcount = 0; nzcount = 0; 1674 for (i=0; i<mbs; i++) { 1675 nmask = 0; 1676 for (j=0; j<bs; j++) { 1677 kmax = rowlengths[rowcount]; 1678 for (k=0; k<kmax; k++) { 1679 tmp = jj[nzcount++]/bs; 1680 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1681 } 1682 rowcount++; 1683 } 1684 browlengths[i] += nmask; 1685 /* zero out the mask elements we set */ 1686 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1687 } 1688 1689 /* create our matrix */ 1690 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1691 B = *A; 1692 a = (Mat_SeqBAIJ*)B->data; 1693 1694 /* set matrix "i" values */ 1695 a->i[0] = 0; 1696 for (i=1; i<= mbs; i++) { 1697 a->i[i] = a->i[i-1] + browlengths[i-1]; 1698 a->ilen[i-1] = browlengths[i-1]; 1699 } 1700 a->nz = 0; 1701 for (i=0; i<mbs; i++) a->nz += browlengths[i]; 1702 1703 /* read in nonzero values */ 1704 ierr = PetscMalloc((nz+extra_rows)*sizeof(Scalar),&aa);CHKERRQ(ierr); 1705 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1706 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1707 1708 /* set "a" and "j" values into matrix */ 1709 nzcount = 0; jcount = 0; 1710 for (i=0; i<mbs; i++) { 1711 nzcountb = nzcount; 1712 nmask = 0; 1713 for (j=0; j<bs; j++) { 1714 kmax = rowlengths[i*bs+j]; 1715 for (k=0; k<kmax; k++) { 1716 tmp = jj[nzcount++]/bs; 1717 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1718 } 1719 rowcount++; 1720 } 1721 /* sort the masked values */ 1722 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1723 1724 /* set "j" values into matrix */ 1725 maskcount = 1; 1726 for (j=0; j<nmask; j++) { 1727 a->j[jcount++] = masked[j]; 1728 mask[masked[j]] = maskcount++; 1729 } 1730 /* set "a" values into matrix */ 1731 ishift = bs2*a->i[i]; 1732 for (j=0; j<bs; j++) { 1733 kmax = rowlengths[i*bs+j]; 1734 for (k=0; k<kmax; k++) { 1735 tmp = jj[nzcountb]/bs ; 1736 block = mask[tmp] - 1; 1737 point = jj[nzcountb] - bs*tmp; 1738 idx = ishift + bs2*block + j + bs*point; 1739 a->a[idx] = (MatScalar)aa[nzcountb++]; 1740 } 1741 } 1742 /* zero out the mask elements we set */ 1743 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1744 } 1745 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1746 1747 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1748 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1749 ierr = PetscFree(aa);CHKERRQ(ierr); 1750 ierr = PetscFree(jj);CHKERRQ(ierr); 1751 ierr = PetscFree(mask);CHKERRQ(ierr); 1752 1753 B->assembled = PETSC_TRUE; 1754 1755 ierr = MatView_Private(B);CHKERRQ(ierr); 1756 PetscFunctionReturn(0); 1757 } 1758 EXTERN_C_END 1759 1760 #undef __FUNC__ 1761 #define __FUNC__ "MatCreateSeqBAIJ" 1762 /*@C 1763 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1764 compressed row) format. For good matrix assembly performance the 1765 user should preallocate the matrix storage by setting the parameter nz 1766 (or the array nnz). By setting these parameters accurately, performance 1767 during matrix assembly can be increased by more than a factor of 50. 1768 1769 Collective on MPI_Comm 1770 1771 Input Parameters: 1772 + comm - MPI communicator, set to PETSC_COMM_SELF 1773 . bs - size of block 1774 . m - number of rows 1775 . n - number of columns 1776 . nz - number of nonzero blocks per block row (same for all rows) 1777 - nnz - array containing the number of nonzero blocks in the various block rows 1778 (possibly different for each block row) or PETSC_NULL 1779 1780 Output Parameter: 1781 . A - the matrix 1782 1783 Options Database Keys: 1784 . -mat_no_unroll - uses code that does not unroll the loops in the 1785 block calculations (much slower) 1786 . -mat_block_size - size of the blocks to use 1787 1788 Level: intermediate 1789 1790 Notes: 1791 A nonzero block is any block that as 1 or more nonzeros in it 1792 1793 The block AIJ format is fully compatible with standard Fortran 77 1794 storage. That is, the stored row and column indices can begin at 1795 either one (as in Fortran) or zero. See the users' manual for details. 1796 1797 Specify the preallocated storage with either nz or nnz (not both). 1798 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1799 allocation. For additional details, see the users manual chapter on 1800 matrices. 1801 1802 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1803 @*/ 1804 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1805 { 1806 int ierr; 1807 1808 PetscFunctionBegin; 1809 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1810 ierr = MatSetType(*A,MATSEQBAIJ);CHKERRQ(ierr); 1811 ierr = MatSeqBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNC__ 1816 #define __FUNC__ "MatSeqBAIJSetPreallocation" 1817 /*@C 1818 MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros 1819 per row in the matrix. For good matrix assembly performance the 1820 user should preallocate the matrix storage by setting the parameter nz 1821 (or the array nnz). By setting these parameters accurately, performance 1822 during matrix assembly can be increased by more than a factor of 50. 1823 1824 Collective on MPI_Comm 1825 1826 Input Parameters: 1827 + A - the matrix 1828 . bs - size of block 1829 . nz - number of block nonzeros per block row (same for all rows) 1830 - nnz - array containing the number of block nonzeros in the various block rows 1831 (possibly different for each block row) or PETSC_NULL 1832 1833 Options Database Keys: 1834 . -mat_no_unroll - uses code that does not unroll the loops in the 1835 block calculations (much slower) 1836 . -mat_block_size - size of the blocks to use 1837 1838 Level: intermediate 1839 1840 Notes: 1841 The block AIJ format is fully compatible with standard Fortran 77 1842 storage. That is, the stored row and column indices can begin at 1843 either one (as in Fortran) or zero. See the users' manual for details. 1844 1845 Specify the preallocated storage with either nz or nnz (not both). 1846 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1847 allocation. For additional details, see the users manual chapter on 1848 matrices. 1849 1850 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1851 @*/ 1852 int MatSeqBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1853 { 1854 Mat_SeqBAIJ *b; 1855 int i,len,ierr,mbs,nbs,bs2,newbs = bs; 1856 PetscTruth flg; 1857 1858 PetscFunctionBegin; 1859 ierr = PetscTypeCompare((PetscObject)B,MATSEQBAIJ,&flg);CHKERRQ(ierr); 1860 if (!flg) PetscFunctionReturn(0); 1861 1862 B->preallocated = PETSC_TRUE; 1863 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&newbs,PETSC_NULL);CHKERRQ(ierr); 1864 if (nnz && newbs != bs) { 1865 SETERRQ(1,"Cannot change blocksize from command line if setting nnz"); 1866 } 1867 bs = newbs; 1868 1869 mbs = B->m/bs; 1870 nbs = B->n/bs; 1871 bs2 = bs*bs; 1872 1873 if (mbs*bs!=B->m || nbs*bs!=B->n) { 1874 SETERRQ3(PETSC_ERR_ARG_SIZ,"Number rows %d, cols %d must be divisible by blocksize %d",B->m,B->n,bs); 1875 } 1876 1877 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 1878 if (nnz) { 1879 for (i=0; i<mbs; i++) { 1880 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1881 } 1882 } 1883 1884 b = (Mat_SeqBAIJ*)B->data; 1885 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1886 if (!flg) { 1887 switch (bs) { 1888 case 1: 1889 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1890 B->ops->solve = MatSolve_SeqBAIJ_1; 1891 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 1892 B->ops->mult = MatMult_SeqBAIJ_1; 1893 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1894 break; 1895 case 2: 1896 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1897 B->ops->solve = MatSolve_SeqBAIJ_2; 1898 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 1899 B->ops->mult = MatMult_SeqBAIJ_2; 1900 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1901 break; 1902 case 3: 1903 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1904 B->ops->solve = MatSolve_SeqBAIJ_3; 1905 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3; 1906 B->ops->mult = MatMult_SeqBAIJ_3; 1907 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1908 break; 1909 case 4: 1910 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1911 B->ops->solve = MatSolve_SeqBAIJ_4; 1912 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_4; 1913 B->ops->mult = MatMult_SeqBAIJ_4; 1914 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1915 break; 1916 case 5: 1917 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1918 B->ops->solve = MatSolve_SeqBAIJ_5; 1919 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5; 1920 B->ops->mult = MatMult_SeqBAIJ_5; 1921 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1922 break; 1923 case 6: 1924 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 1925 B->ops->solve = MatSolve_SeqBAIJ_6; 1926 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6; 1927 B->ops->mult = MatMult_SeqBAIJ_6; 1928 B->ops->multadd = MatMultAdd_SeqBAIJ_6; 1929 break; 1930 case 7: 1931 B->ops->mult = MatMult_SeqBAIJ_7; 1932 B->ops->solve = MatSolve_SeqBAIJ_7; 1933 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7; 1934 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1935 break; 1936 default: 1937 B->ops->mult = MatMult_SeqBAIJ_N; 1938 B->ops->solve = MatSolve_SeqBAIJ_N; 1939 B->ops->multadd = MatMultAdd_SeqBAIJ_N; 1940 break; 1941 } 1942 } 1943 b->bs = bs; 1944 b->mbs = mbs; 1945 b->nbs = nbs; 1946 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 1947 if (!nnz) { 1948 if (nz == PETSC_DEFAULT) nz = 5; 1949 else if (nz <= 0) nz = 1; 1950 for (i=0; i<mbs; i++) b->imax[i] = nz; 1951 nz = nz*mbs; 1952 } else { 1953 nz = 0; 1954 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1955 } 1956 1957 /* allocate the matrix space */ 1958 len = nz*sizeof(int) + nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1959 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 1960 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1961 b->j = (int*)(b->a + nz*bs2); 1962 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 1963 b->i = b->j + nz; 1964 b->singlemalloc = PETSC_TRUE; 1965 1966 b->i[0] = 0; 1967 for (i=1; i<mbs+1; i++) { 1968 b->i[i] = b->i[i-1] + b->imax[i-1]; 1969 } 1970 1971 /* b->ilen will count nonzeros in each block row so far. */ 1972 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->ilen); 1973 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1974 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1975 1976 b->bs = bs; 1977 b->bs2 = bs2; 1978 b->mbs = mbs; 1979 b->nz = 0; 1980 b->maxnz = nz*bs2; 1981 B->info.nz_unneeded = (PetscReal)b->maxnz; 1982 PetscFunctionReturn(0); 1983 } 1984 1985