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