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