1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the SBAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "../src/mat/impls/baij/seq/baij.h" /*I "petscmat.h" I*/ 8 #include "../src/mat/impls/sbaij/seq/sbaij.h" 9 10 #include "../src/mat/impls/sbaij/seq/relax.h" 11 #define USESHORT 12 #include "../src/mat/impls/sbaij/seq/relax.h" 13 14 extern PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat,PetscTruth); 15 #define CHUNKSIZE 10 16 17 18 /* 19 Checks for missing diagonals 20 */ 21 #undef __FUNCT__ 22 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 23 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscTruth *missing,PetscInt *dd) 24 { 25 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 26 PetscErrorCode ierr; 27 PetscInt *diag,*jj = a->j,i; 28 29 PetscFunctionBegin; 30 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 31 diag = a->diag; 32 *missing = PETSC_FALSE; 33 for (i=0; i<a->mbs; i++) { 34 if (jj[diag[i]] != i) { 35 *missing = PETSC_TRUE; 36 if (dd) *dd = i; 37 break; 38 } 39 } 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 45 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 46 { 47 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 48 PetscErrorCode ierr; 49 PetscInt i; 50 51 PetscFunctionBegin; 52 if (!a->diag) { 53 ierr = PetscMalloc(a->mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 54 ierr = PetscLogObjectMemory(A,a->mbs*sizeof(PetscInt));CHKERRQ(ierr); 55 a->free_diag = PETSC_TRUE; 56 } 57 for (i=0; i<a->mbs; i++) a->diag[i] = a->i[i]; 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 63 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 64 { 65 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 66 PetscInt i,j,n = a->mbs,nz = a->i[n],bs = A->rmap->bs; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 *nn = n; 71 if (!ia) PetscFunctionReturn(0); 72 if (!blockcompressed) { 73 /* malloc & create the natural set of indices */ 74 ierr = PetscMalloc2((n+1)*bs,PetscInt,ia,nz*bs,PetscInt,ja);CHKERRQ(ierr); 75 for (i=0; i<n+1; i++) { 76 for (j=0; j<bs; j++) { 77 *ia[i*bs+j] = a->i[i]*bs+j+oshift; 78 } 79 } 80 for (i=0; i<nz; i++) { 81 for (j=0; j<bs; j++) { 82 *ja[i*bs+j] = a->j[i]*bs+j+oshift; 83 } 84 } 85 } else { /* blockcompressed */ 86 if (oshift == 1) { 87 /* temporarily add 1 to i and j indices */ 88 for (i=0; i<nz; i++) a->j[i]++; 89 for (i=0; i<n+1; i++) a->i[i]++; 90 } 91 *ia = a->i; *ja = a->j; 92 } 93 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 99 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 100 { 101 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 102 PetscInt i,n = a->mbs,nz = a->i[n]; 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 if (!ia) PetscFunctionReturn(0); 107 108 if (!blockcompressed) { 109 ierr = PetscFree2(*ia,*ja);CHKERRQ(ierr); 110 } else if (oshift == 1) { /* blockcompressed */ 111 for (i=0; i<nz; i++) a->j[i]--; 112 for (i=0; i<n+1; i++) a->i[i]--; 113 } 114 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 120 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 121 { 122 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 #if defined(PETSC_USE_LOG) 127 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz); 128 #endif 129 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 130 if (a->free_diag){ierr = PetscFree(a->diag);CHKERRQ(ierr);} 131 if (a->row) {ierr = ISDestroy(a->row);CHKERRQ(ierr);} 132 if (a->col){ierr = ISDestroy(a->col);CHKERRQ(ierr);} 133 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 134 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 135 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 136 if (a->free_diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 137 if (a->free_imax_ilen) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);} 138 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 139 ierr = PetscFree(a->sor_work);CHKERRQ(ierr); 140 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 141 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 142 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 143 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 144 if (a->free_jshort) {ierr = PetscFree(a->jshort);CHKERRQ(ierr);} 145 ierr = PetscFree(a->inew);CHKERRQ(ierr); 146 if (a->parent) {ierr = MatDestroy(a->parent);CHKERRQ(ierr);} 147 ierr = PetscFree(a);CHKERRQ(ierr); 148 149 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 150 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 151 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 152 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 153 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 154 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 155 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 161 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscTruth flg) 162 { 163 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 switch (op) { 168 case MAT_ROW_ORIENTED: 169 a->roworiented = flg; 170 break; 171 case MAT_KEEP_NONZERO_PATTERN: 172 a->keepnonzeropattern = flg; 173 break; 174 case MAT_NEW_NONZERO_LOCATIONS: 175 a->nonew = (flg ? 0 : 1); 176 break; 177 case MAT_NEW_NONZERO_LOCATION_ERR: 178 a->nonew = (flg ? -1 : 0); 179 break; 180 case MAT_NEW_NONZERO_ALLOCATION_ERR: 181 a->nonew = (flg ? -2 : 0); 182 break; 183 case MAT_UNUSED_NONZERO_LOCATION_ERR: 184 a->nounused = (flg ? -1 : 0); 185 break; 186 case MAT_NEW_DIAGONALS: 187 case MAT_IGNORE_OFF_PROC_ENTRIES: 188 case MAT_USE_HASH_TABLE: 189 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 190 break; 191 case MAT_HERMITIAN: 192 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call MatAssemblyEnd() first"); 193 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 194 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort; 195 } else if (A->cmap->bs == 1) { 196 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian; 197 } else SETERRQ(PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1"); 198 break; 199 case MAT_SYMMETRIC: 200 case MAT_STRUCTURALLY_SYMMETRIC: 201 case MAT_SYMMETRY_ETERNAL: 202 if (!flg) SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric"); 203 ierr = PetscInfo1(A,"Option %s not relevent\n",MatOptions[op]);CHKERRQ(ierr); 204 break; 205 case MAT_IGNORE_LOWER_TRIANGULAR: 206 a->ignore_ltriangular = flg; 207 break; 208 case MAT_ERROR_LOWER_TRIANGULAR: 209 a->ignore_ltriangular = flg; 210 break; 211 case MAT_GETROW_UPPERTRIANGULAR: 212 a->getrow_utriangular = flg; 213 break; 214 default: 215 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 216 } 217 PetscFunctionReturn(0); 218 } 219 220 #undef __FUNCT__ 221 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 222 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v) 223 { 224 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 225 PetscErrorCode ierr; 226 PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 227 MatScalar *aa,*aa_i; 228 PetscScalar *v_i; 229 230 PetscFunctionBegin; 231 if (A && !a->getrow_utriangular) SETERRQ(PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()"); 232 /* Get the upper triangular part of the row */ 233 bs = A->rmap->bs; 234 ai = a->i; 235 aj = a->j; 236 aa = a->a; 237 bs2 = a->bs2; 238 239 if (row < 0 || row >= A->rmap->N) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row); 240 241 bn = row/bs; /* Block number */ 242 bp = row % bs; /* Block position */ 243 M = ai[bn+1] - ai[bn]; 244 *ncols = bs*M; 245 246 if (v) { 247 *v = 0; 248 if (*ncols) { 249 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 250 for (i=0; i<M; i++) { /* for each block in the block row */ 251 v_i = *v + i*bs; 252 aa_i = aa + bs2*(ai[bn] + i); 253 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 254 } 255 } 256 } 257 258 if (cols) { 259 *cols = 0; 260 if (*ncols) { 261 ierr = PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);CHKERRQ(ierr); 262 for (i=0; i<M; i++) { /* for each block in the block row */ 263 cols_i = *cols + i*bs; 264 itmp = bs*aj[ai[bn] + i]; 265 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 266 } 267 } 268 } 269 270 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 271 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 272 #ifdef column_search 273 v_i = *v + M*bs; 274 cols_i = *cols + M*bs; 275 for (i=0; i<bn; i++){ /* for each block row */ 276 M = ai[i+1] - ai[i]; 277 for (j=0; j<M; j++){ 278 itmp = aj[ai[i] + j]; /* block column value */ 279 if (itmp == bn){ 280 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 281 for (k=0; k<bs; k++) { 282 *cols_i++ = i*bs+k; 283 *v_i++ = aa_i[k]; 284 } 285 *ncols += bs; 286 break; 287 } 288 } 289 } 290 #endif 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 296 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 297 { 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 302 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 303 PetscFunctionReturn(0); 304 } 305 306 #undef __FUNCT__ 307 #define __FUNCT__ "MatGetRowUpperTriangular_SeqSBAIJ" 308 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 309 { 310 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 311 312 PetscFunctionBegin; 313 a->getrow_utriangular = PETSC_TRUE; 314 PetscFunctionReturn(0); 315 } 316 #undef __FUNCT__ 317 #define __FUNCT__ "MatRestoreRowUpperTriangular_SeqSBAIJ" 318 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 319 { 320 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 321 322 PetscFunctionBegin; 323 a->getrow_utriangular = PETSC_FALSE; 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 329 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 330 { 331 PetscErrorCode ierr; 332 PetscFunctionBegin; 333 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 334 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 335 } 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 341 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 342 { 343 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 344 PetscErrorCode ierr; 345 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 346 const char *name; 347 PetscViewerFormat format; 348 349 PetscFunctionBegin; 350 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 351 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 352 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 353 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 354 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 355 Mat aij; 356 357 if (A->factor && bs>1){ 358 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 362 ierr = MatView(aij,viewer);CHKERRQ(ierr); 363 ierr = MatDestroy(aij);CHKERRQ(ierr); 364 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 365 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 366 for (i=0; i<a->mbs; i++) { 367 for (j=0; j<bs; j++) { 368 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 369 for (k=a->i[i]; k<a->i[i+1]; k++) { 370 for (l=0; l<bs; l++) { 371 #if defined(PETSC_USE_COMPLEX) 372 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 373 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 374 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 375 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 376 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 377 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 378 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 379 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 380 } 381 #else 382 if (a->a[bs2*k + l*bs + j] != 0.0) { 383 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 384 } 385 #endif 386 } 387 } 388 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 389 } 390 } 391 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 392 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 393 PetscFunctionReturn(0); 394 } else { 395 if (A->factor && bs>1){ 396 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored. MatView_SeqSBAIJ_ASCII() may not display complete or logically correct entries!\n");CHKERRQ(ierr); 397 } 398 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 399 for (i=0; i<a->mbs; i++) { 400 for (j=0; j<bs; j++) { 401 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 402 for (k=a->i[i]; k<a->i[i+1]; k++) { 403 for (l=0; l<bs; l++) { 404 #if defined(PETSC_USE_COMPLEX) 405 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 406 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 407 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 408 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 409 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 410 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 411 } else { 412 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 413 } 414 #else 415 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 416 #endif 417 } 418 } 419 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 420 } 421 } 422 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 423 } 424 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 425 PetscFunctionReturn(0); 426 } 427 428 #undef __FUNCT__ 429 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 430 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 431 { 432 Mat A = (Mat) Aa; 433 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 434 PetscErrorCode ierr; 435 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 436 PetscMPIInt rank; 437 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 438 MatScalar *aa; 439 MPI_Comm comm; 440 PetscViewer viewer; 441 442 PetscFunctionBegin; 443 /* 444 This is nasty. If this is called from an originally parallel matrix 445 then all processes call this,but only the first has the matrix so the 446 rest should return immediately. 447 */ 448 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 449 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 450 if (rank) PetscFunctionReturn(0); 451 452 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 453 454 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 455 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 456 457 /* loop over matrix elements drawing boxes */ 458 color = PETSC_DRAW_BLUE; 459 for (i=0,row=0; i<mbs; i++,row+=bs) { 460 for (j=a->i[i]; j<a->i[i+1]; j++) { 461 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 462 x_l = a->j[j]*bs; x_r = x_l + 1.0; 463 aa = a->a + j*bs2; 464 for (k=0; k<bs; k++) { 465 for (l=0; l<bs; l++) { 466 if (PetscRealPart(*aa++) >= 0.) continue; 467 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 468 } 469 } 470 } 471 } 472 color = PETSC_DRAW_CYAN; 473 for (i=0,row=0; i<mbs; i++,row+=bs) { 474 for (j=a->i[i]; j<a->i[i+1]; j++) { 475 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 476 x_l = a->j[j]*bs; x_r = x_l + 1.0; 477 aa = a->a + j*bs2; 478 for (k=0; k<bs; k++) { 479 for (l=0; l<bs; l++) { 480 if (PetscRealPart(*aa++) != 0.) continue; 481 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 482 } 483 } 484 } 485 } 486 487 color = PETSC_DRAW_RED; 488 for (i=0,row=0; i<mbs; i++,row+=bs) { 489 for (j=a->i[i]; j<a->i[i+1]; j++) { 490 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 491 x_l = a->j[j]*bs; x_r = x_l + 1.0; 492 aa = a->a + j*bs2; 493 for (k=0; k<bs; k++) { 494 for (l=0; l<bs; l++) { 495 if (PetscRealPart(*aa++) <= 0.) continue; 496 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 497 } 498 } 499 } 500 } 501 PetscFunctionReturn(0); 502 } 503 504 #undef __FUNCT__ 505 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 506 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 507 { 508 PetscErrorCode ierr; 509 PetscReal xl,yl,xr,yr,w,h; 510 PetscDraw draw; 511 PetscTruth isnull; 512 513 PetscFunctionBegin; 514 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 515 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 516 517 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 518 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 519 xr += w; yr += h; xl = -w; yl = -h; 520 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 521 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 522 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 523 PetscFunctionReturn(0); 524 } 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "MatView_SeqSBAIJ" 528 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 529 { 530 PetscErrorCode ierr; 531 PetscTruth iascii,isdraw; 532 FILE *file = 0; 533 534 PetscFunctionBegin; 535 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 536 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 537 if (iascii){ 538 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 539 } else if (isdraw) { 540 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 541 } else { 542 Mat B; 543 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 544 ierr = MatView(B,viewer);CHKERRQ(ierr); 545 ierr = MatDestroy(B);CHKERRQ(ierr); 546 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 547 if (file) { 548 fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs); 549 } 550 } 551 PetscFunctionReturn(0); 552 } 553 554 555 #undef __FUNCT__ 556 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 557 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 558 { 559 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 560 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 561 PetscInt *ai = a->i,*ailen = a->ilen; 562 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 563 MatScalar *ap,*aa = a->a; 564 565 PetscFunctionBegin; 566 for (k=0; k<m; k++) { /* loop over rows */ 567 row = im[k]; brow = row/bs; 568 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 569 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 570 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 571 nrow = ailen[brow]; 572 for (l=0; l<n; l++) { /* loop over columns */ 573 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 574 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 575 col = in[l] ; 576 bcol = col/bs; 577 cidx = col%bs; 578 ridx = row%bs; 579 high = nrow; 580 low = 0; /* assume unsorted */ 581 while (high-low > 5) { 582 t = (low+high)/2; 583 if (rp[t] > bcol) high = t; 584 else low = t; 585 } 586 for (i=low; i<high; i++) { 587 if (rp[i] > bcol) break; 588 if (rp[i] == bcol) { 589 *v++ = ap[bs2*i+bs*cidx+ridx]; 590 goto finished; 591 } 592 } 593 *v++ = 0.0; 594 finished:; 595 } 596 } 597 PetscFunctionReturn(0); 598 } 599 600 601 #undef __FUNCT__ 602 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 603 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 604 { 605 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 606 PetscErrorCode ierr; 607 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 608 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 609 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 610 PetscTruth roworiented=a->roworiented; 611 const PetscScalar *value = v; 612 MatScalar *ap,*aa = a->a,*bap; 613 614 PetscFunctionBegin; 615 if (roworiented) { 616 stepval = (n-1)*bs; 617 } else { 618 stepval = (m-1)*bs; 619 } 620 for (k=0; k<m; k++) { /* loop over added rows */ 621 row = im[k]; 622 if (row < 0) continue; 623 #if defined(PETSC_USE_DEBUG) 624 if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 625 #endif 626 rp = aj + ai[row]; 627 ap = aa + bs2*ai[row]; 628 rmax = imax[row]; 629 nrow = ailen[row]; 630 low = 0; 631 high = nrow; 632 for (l=0; l<n; l++) { /* loop over added columns */ 633 if (in[l] < 0) continue; 634 col = in[l]; 635 #if defined(PETSC_USE_DEBUG) 636 if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 637 #endif 638 if (col < row) { 639 if (a->ignore_ltriangular) { 640 continue; /* ignore lower triangular block */ 641 } else { 642 SETERRQ(PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 643 } 644 } 645 if (roworiented) { 646 value = v + k*(stepval+bs)*bs + l*bs; 647 } else { 648 value = v + l*(stepval+bs)*bs + k*bs; 649 } 650 if (col <= lastcol) low = 0; else high = nrow; 651 lastcol = col; 652 while (high-low > 7) { 653 t = (low+high)/2; 654 if (rp[t] > col) high = t; 655 else low = t; 656 } 657 for (i=low; i<high; i++) { 658 if (rp[i] > col) break; 659 if (rp[i] == col) { 660 bap = ap + bs2*i; 661 if (roworiented) { 662 if (is == ADD_VALUES) { 663 for (ii=0; ii<bs; ii++,value+=stepval) { 664 for (jj=ii; jj<bs2; jj+=bs) { 665 bap[jj] += *value++; 666 } 667 } 668 } else { 669 for (ii=0; ii<bs; ii++,value+=stepval) { 670 for (jj=ii; jj<bs2; jj+=bs) { 671 bap[jj] = *value++; 672 } 673 } 674 } 675 } else { 676 if (is == ADD_VALUES) { 677 for (ii=0; ii<bs; ii++,value+=stepval) { 678 for (jj=0; jj<bs; jj++) { 679 *bap++ += *value++; 680 } 681 } 682 } else { 683 for (ii=0; ii<bs; ii++,value+=stepval) { 684 for (jj=0; jj<bs; jj++) { 685 *bap++ = *value++; 686 } 687 } 688 } 689 } 690 goto noinsert2; 691 } 692 } 693 if (nonew == 1) goto noinsert2; 694 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 695 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 696 N = nrow++ - 1; high++; 697 /* shift up all the later entries in this row */ 698 for (ii=N; ii>=i; ii--) { 699 rp[ii+1] = rp[ii]; 700 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 701 } 702 if (N >= i) { 703 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 704 } 705 rp[i] = col; 706 bap = ap + bs2*i; 707 if (roworiented) { 708 for (ii=0; ii<bs; ii++,value+=stepval) { 709 for (jj=ii; jj<bs2; jj+=bs) { 710 bap[jj] = *value++; 711 } 712 } 713 } else { 714 for (ii=0; ii<bs; ii++,value+=stepval) { 715 for (jj=0; jj<bs; jj++) { 716 *bap++ = *value++; 717 } 718 } 719 } 720 noinsert2:; 721 low = i; 722 } 723 ailen[row] = nrow; 724 } 725 PetscFunctionReturn(0); 726 } 727 728 /* 729 This is not yet used 730 */ 731 #undef __FUNCT__ 732 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 733 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 734 { 735 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 736 PetscErrorCode ierr; 737 const PetscInt *ai = a->i, *aj = a->j,*cols; 738 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 739 PetscTruth flag; 740 741 PetscFunctionBegin; 742 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 743 while (i < m){ 744 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 745 /* Limits the number of elements in a node to 'a->inode.limit' */ 746 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 747 nzy = ai[j+1] - ai[j]; 748 if (nzy != (nzx - j + i)) break; 749 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 750 if (!flag) break; 751 } 752 ns[node_count++] = blk_size; 753 i = j; 754 } 755 if (!a->inode.size && m && node_count > .9*m) { 756 ierr = PetscFree(ns);CHKERRQ(ierr); 757 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 758 } else { 759 a->inode.node_count = node_count; 760 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 761 ierr = PetscLogObjectMemory(A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 762 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 763 ierr = PetscFree(ns);CHKERRQ(ierr); 764 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 765 766 /* count collections of adjacent columns in each inode */ 767 row = 0; 768 cnt = 0; 769 for (i=0; i<node_count; i++) { 770 cols = aj + ai[row] + a->inode.size[i]; 771 nz = ai[row+1] - ai[row] - a->inode.size[i]; 772 for (j=1; j<nz; j++) { 773 if (cols[j] != cols[j-1]+1) { 774 cnt++; 775 } 776 } 777 cnt++; 778 row += a->inode.size[i]; 779 } 780 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 781 cnt = 0; 782 row = 0; 783 for (i=0; i<node_count; i++) { 784 cols = aj + ai[row] + a->inode.size[i]; 785 CHKMEMQ; 786 counts[2*cnt] = cols[0]; 787 CHKMEMQ; 788 nz = ai[row+1] - ai[row] - a->inode.size[i]; 789 cnt2 = 1; 790 for (j=1; j<nz; j++) { 791 if (cols[j] != cols[j-1]+1) { 792 CHKMEMQ; 793 counts[2*(cnt++)+1] = cnt2; 794 counts[2*cnt] = cols[j]; 795 CHKMEMQ; 796 cnt2 = 1; 797 } else cnt2++; 798 } 799 CHKMEMQ; 800 counts[2*(cnt++)+1] = cnt2; 801 CHKMEMQ; 802 row += a->inode.size[i]; 803 } 804 ierr = PetscIntView(2*cnt,counts,0); 805 } 806 PetscFunctionReturn(0); 807 } 808 809 #undef __FUNCT__ 810 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 811 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 812 { 813 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 814 PetscErrorCode ierr; 815 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 816 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 817 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 818 MatScalar *aa = a->a,*ap; 819 820 PetscFunctionBegin; 821 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 822 823 if (m) rmax = ailen[0]; 824 for (i=1; i<mbs; i++) { 825 /* move each row back by the amount of empty slots (fshift) before it*/ 826 fshift += imax[i-1] - ailen[i-1]; 827 rmax = PetscMax(rmax,ailen[i]); 828 if (fshift) { 829 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 830 N = ailen[i]; 831 for (j=0; j<N; j++) { 832 ip[j-fshift] = ip[j]; 833 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 834 } 835 } 836 ai[i] = ai[i-1] + ailen[i-1]; 837 } 838 if (mbs) { 839 fshift += imax[mbs-1] - ailen[mbs-1]; 840 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 841 } 842 /* reset ilen and imax for each row */ 843 for (i=0; i<mbs; i++) { 844 ailen[i] = imax[i] = ai[i+1] - ai[i]; 845 } 846 a->nz = ai[mbs]; 847 848 /* diagonals may have moved, reset it */ 849 if (a->diag) { 850 ierr = PetscMemcpy(a->diag,ai,mbs*sizeof(PetscInt));CHKERRQ(ierr); 851 } 852 if (fshift && a->nounused == -1) { 853 SETERRQ4(PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2); 854 } 855 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr); 856 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 857 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 858 a->reallocs = 0; 859 A->info.nz_unneeded = (PetscReal)fshift*bs2; 860 a->idiagvalid = PETSC_FALSE; 861 862 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 863 if (!a->jshort) { 864 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 865 ierr = PetscLogObjectMemory(A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 866 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 867 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 868 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 869 a->free_jshort = PETSC_TRUE; 870 } 871 } 872 PetscFunctionReturn(0); 873 } 874 875 /* 876 This function returns an array of flags which indicate the locations of contiguous 877 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 878 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 879 Assume: sizes should be long enough to hold all the values. 880 */ 881 #undef __FUNCT__ 882 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 883 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 884 { 885 PetscInt i,j,k,row; 886 PetscTruth flg; 887 888 PetscFunctionBegin; 889 for (i=0,j=0; i<n; j++) { 890 row = idx[i]; 891 if (row%bs!=0) { /* Not the begining of a block */ 892 sizes[j] = 1; 893 i++; 894 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 895 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 896 i++; 897 } else { /* Begining of the block, so check if the complete block exists */ 898 flg = PETSC_TRUE; 899 for (k=1; k<bs; k++) { 900 if (row+k != idx[i+k]) { /* break in the block */ 901 flg = PETSC_FALSE; 902 break; 903 } 904 } 905 if (flg) { /* No break in the bs */ 906 sizes[j] = bs; 907 i+= bs; 908 } else { 909 sizes[j] = 1; 910 i++; 911 } 912 } 913 } 914 *bs_max = j; 915 PetscFunctionReturn(0); 916 } 917 918 919 /* Only add/insert a(i,j) with i<=j (blocks). 920 Any a(i,j) with i>j input by user is ingored. 921 */ 922 923 #undef __FUNCT__ 924 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 925 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 926 { 927 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 928 PetscErrorCode ierr; 929 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 930 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 931 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 932 PetscInt ridx,cidx,bs2=a->bs2; 933 MatScalar *ap,value,*aa=a->a,*bap; 934 935 PetscFunctionBegin; 936 if (v) PetscValidScalarPointer(v,6); 937 for (k=0; k<m; k++) { /* loop over added rows */ 938 row = im[k]; /* row number */ 939 brow = row/bs; /* block row number */ 940 if (row < 0) continue; 941 #if defined(PETSC_USE_DEBUG) 942 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 943 #endif 944 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 945 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 946 rmax = imax[brow]; /* maximum space allocated for this row */ 947 nrow = ailen[brow]; /* actual length of this row */ 948 low = 0; 949 950 for (l=0; l<n; l++) { /* loop over added columns */ 951 if (in[l] < 0) continue; 952 #if defined(PETSC_USE_DEBUG) 953 if (in[l] >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->N-1); 954 #endif 955 col = in[l]; 956 bcol = col/bs; /* block col number */ 957 958 if (brow > bcol) { 959 if (a->ignore_ltriangular){ 960 continue; /* ignore lower triangular values */ 961 } else { 962 SETERRQ(PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 963 } 964 } 965 966 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 967 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 968 /* element value a(k,l) */ 969 if (roworiented) { 970 value = v[l + k*n]; 971 } else { 972 value = v[k + l*m]; 973 } 974 975 /* move pointer bap to a(k,l) quickly and add/insert value */ 976 if (col <= lastcol) low = 0; high = nrow; 977 lastcol = col; 978 while (high-low > 7) { 979 t = (low+high)/2; 980 if (rp[t] > bcol) high = t; 981 else low = t; 982 } 983 for (i=low; i<high; i++) { 984 if (rp[i] > bcol) break; 985 if (rp[i] == bcol) { 986 bap = ap + bs2*i + bs*cidx + ridx; 987 if (is == ADD_VALUES) *bap += value; 988 else *bap = value; 989 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 990 if (brow == bcol && ridx < cidx){ 991 bap = ap + bs2*i + bs*ridx + cidx; 992 if (is == ADD_VALUES) *bap += value; 993 else *bap = value; 994 } 995 goto noinsert1; 996 } 997 } 998 999 if (nonew == 1) goto noinsert1; 1000 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1001 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1002 1003 N = nrow++ - 1; high++; 1004 /* shift up all the later entries in this row */ 1005 for (ii=N; ii>=i; ii--) { 1006 rp[ii+1] = rp[ii]; 1007 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1008 } 1009 if (N>=i) { 1010 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1011 } 1012 rp[i] = bcol; 1013 ap[bs2*i + bs*cidx + ridx] = value; 1014 noinsert1:; 1015 low = i; 1016 } 1017 } /* end of loop over added columns */ 1018 ailen[brow] = nrow; 1019 } /* end of loop over added rows */ 1020 PetscFunctionReturn(0); 1021 } 1022 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 1025 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1026 { 1027 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1028 Mat outA; 1029 PetscErrorCode ierr; 1030 PetscTruth row_identity; 1031 1032 PetscFunctionBegin; 1033 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1034 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1035 if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported"); 1036 if (inA->rmap->bs != 1) SETERRQ1(PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 1037 1038 outA = inA; 1039 inA->factor = MAT_FACTOR_ICC; 1040 1041 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1042 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1043 1044 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1045 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr); } 1046 a->row = row; 1047 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1048 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr); } 1049 a->col = row; 1050 1051 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1052 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1053 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1054 1055 if (!a->solve_work) { 1056 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1057 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1058 } 1059 1060 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1061 PetscFunctionReturn(0); 1062 } 1063 1064 EXTERN_C_BEGIN 1065 #undef __FUNCT__ 1066 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1067 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1068 { 1069 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1070 PetscInt i,nz,n; 1071 1072 PetscFunctionBegin; 1073 nz = baij->maxnz; 1074 n = mat->cmap->n; 1075 for (i=0; i<nz; i++) { 1076 baij->j[i] = indices[i]; 1077 } 1078 baij->nz = nz; 1079 for (i=0; i<n; i++) { 1080 baij->ilen[i] = baij->imax[i]; 1081 } 1082 PetscFunctionReturn(0); 1083 } 1084 EXTERN_C_END 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1088 /*@ 1089 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1090 in the matrix. 1091 1092 Input Parameters: 1093 + mat - the SeqSBAIJ matrix 1094 - indices - the column indices 1095 1096 Level: advanced 1097 1098 Notes: 1099 This can be called if you have precomputed the nonzero structure of the 1100 matrix and want to provide it to the matrix object to improve the performance 1101 of the MatSetValues() operation. 1102 1103 You MUST have set the correct numbers of nonzeros per row in the call to 1104 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1105 1106 MUST be called before any calls to MatSetValues() 1107 1108 .seealso: MatCreateSeqSBAIJ 1109 @*/ 1110 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1111 { 1112 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1113 1114 PetscFunctionBegin; 1115 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1116 PetscValidPointer(indices,2); 1117 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1118 if (f) { 1119 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1120 } else { 1121 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1122 } 1123 PetscFunctionReturn(0); 1124 } 1125 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1128 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1129 { 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 /* If the two matrices have the same copy implementation, use fast copy. */ 1134 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1135 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1136 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1137 1138 if (a->i[A->rmap->N] != b->i[B->rmap->N]) { 1139 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1140 } 1141 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1142 } else { 1143 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1144 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1145 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1146 } 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1152 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1153 { 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1158 PetscFunctionReturn(0); 1159 } 1160 1161 #undef __FUNCT__ 1162 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1163 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1164 { 1165 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1166 PetscFunctionBegin; 1167 *array = a->a; 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1173 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1174 { 1175 PetscFunctionBegin; 1176 PetscFunctionReturn(0); 1177 } 1178 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1181 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1182 { 1183 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1184 PetscErrorCode ierr; 1185 PetscInt i,bs=Y->rmap->bs,bs2,j; 1186 PetscBLASInt one = 1,bnz = PetscBLASIntCast(x->nz); 1187 1188 PetscFunctionBegin; 1189 if (str == SAME_NONZERO_PATTERN) { 1190 PetscScalar alpha = a; 1191 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1192 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1193 if (y->xtoy && y->XtoY != X) { 1194 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1195 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1196 } 1197 if (!y->xtoy) { /* get xtoy */ 1198 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1199 y->XtoY = X; 1200 } 1201 bs2 = bs*bs; 1202 for (i=0; i<x->nz; i++) { 1203 j = 0; 1204 while (j < bs2){ 1205 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1206 j++; 1207 } 1208 } 1209 ierr = PetscInfo3(Y,"ratio of nnz_s(X)/nnz_s(Y): %D/%D = %G\n",bs2*x->nz,bs2*y->nz,(PetscReal)(bs2*x->nz)/(bs2*y->nz));CHKERRQ(ierr); 1210 } else { 1211 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1212 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1213 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1214 } 1215 PetscFunctionReturn(0); 1216 } 1217 1218 #undef __FUNCT__ 1219 #define __FUNCT__ "MatSetBlockSize_SeqSBAIJ" 1220 PetscErrorCode MatSetBlockSize_SeqSBAIJ(Mat A,PetscInt bs) 1221 { 1222 PetscInt rbs,cbs; 1223 PetscErrorCode ierr; 1224 1225 PetscFunctionBegin; 1226 ierr = PetscLayoutGetBlockSize(A->rmap,&rbs);CHKERRQ(ierr); 1227 ierr = PetscLayoutGetBlockSize(A->cmap,&cbs);CHKERRQ(ierr); 1228 if (rbs != bs) SETERRQ2(PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,rbs); 1229 if (cbs != bs) SETERRQ2(PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,cbs); 1230 PetscFunctionReturn(0); 1231 } 1232 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1235 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1236 { 1237 PetscFunctionBegin; 1238 *flg = PETSC_TRUE; 1239 PetscFunctionReturn(0); 1240 } 1241 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1244 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1245 { 1246 PetscFunctionBegin; 1247 *flg = PETSC_TRUE; 1248 PetscFunctionReturn(0); 1249 } 1250 1251 #undef __FUNCT__ 1252 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1253 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1254 { 1255 PetscFunctionBegin; 1256 *flg = PETSC_FALSE; 1257 PetscFunctionReturn(0); 1258 } 1259 1260 #undef __FUNCT__ 1261 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1262 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1263 { 1264 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1265 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1266 MatScalar *aa = a->a; 1267 1268 PetscFunctionBegin; 1269 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1275 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1276 { 1277 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1278 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1279 MatScalar *aa = a->a; 1280 1281 PetscFunctionBegin; 1282 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 /* -------------------------------------------------------------------*/ 1287 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1288 MatGetRow_SeqSBAIJ, 1289 MatRestoreRow_SeqSBAIJ, 1290 MatMult_SeqSBAIJ_N, 1291 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1292 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1293 MatMultAdd_SeqSBAIJ_N, 1294 0, 1295 0, 1296 0, 1297 /*10*/ 0, 1298 0, 1299 MatCholeskyFactor_SeqSBAIJ, 1300 MatSOR_SeqSBAIJ, 1301 MatTranspose_SeqSBAIJ, 1302 /*15*/ MatGetInfo_SeqSBAIJ, 1303 MatEqual_SeqSBAIJ, 1304 MatGetDiagonal_SeqSBAIJ, 1305 MatDiagonalScale_SeqSBAIJ, 1306 MatNorm_SeqSBAIJ, 1307 /*20*/ 0, 1308 MatAssemblyEnd_SeqSBAIJ, 1309 MatSetOption_SeqSBAIJ, 1310 MatZeroEntries_SeqSBAIJ, 1311 /*24*/ 0, 1312 0, 1313 0, 1314 0, 1315 0, 1316 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1317 0, 1318 0, 1319 MatGetArray_SeqSBAIJ, 1320 MatRestoreArray_SeqSBAIJ, 1321 /*34*/ MatDuplicate_SeqSBAIJ, 1322 0, 1323 0, 1324 0, 1325 MatICCFactor_SeqSBAIJ, 1326 /*39*/ MatAXPY_SeqSBAIJ, 1327 MatGetSubMatrices_SeqSBAIJ, 1328 MatIncreaseOverlap_SeqSBAIJ, 1329 MatGetValues_SeqSBAIJ, 1330 MatCopy_SeqSBAIJ, 1331 /*44*/ 0, 1332 MatScale_SeqSBAIJ, 1333 0, 1334 0, 1335 0, 1336 /*49*/ MatSetBlockSize_SeqSBAIJ, 1337 MatGetRowIJ_SeqSBAIJ, 1338 MatRestoreRowIJ_SeqSBAIJ, 1339 0, 1340 0, 1341 /*54*/ 0, 1342 0, 1343 0, 1344 0, 1345 MatSetValuesBlocked_SeqSBAIJ, 1346 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1347 0, 1348 0, 1349 0, 1350 0, 1351 /*64*/ 0, 1352 0, 1353 0, 1354 0, 1355 0, 1356 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1357 0, 1358 0, 1359 0, 1360 0, 1361 /*74*/ 0, 1362 0, 1363 0, 1364 0, 1365 0, 1366 /*79*/ 0, 1367 0, 1368 0, 1369 #if !defined(PETSC_USE_COMPLEX) 1370 MatGetInertia_SeqSBAIJ, 1371 #else 1372 0, 1373 #endif 1374 MatLoad_SeqSBAIJ, 1375 /*84*/ MatIsSymmetric_SeqSBAIJ, 1376 MatIsHermitian_SeqSBAIJ, 1377 MatIsStructurallySymmetric_SeqSBAIJ, 1378 0, 1379 0, 1380 /*89*/ 0, 1381 0, 1382 0, 1383 0, 1384 0, 1385 /*94*/ 0, 1386 0, 1387 0, 1388 0, 1389 0, 1390 /*99*/ 0, 1391 0, 1392 0, 1393 0, 1394 0, 1395 /*104*/0, 1396 MatRealPart_SeqSBAIJ, 1397 MatImaginaryPart_SeqSBAIJ, 1398 MatGetRowUpperTriangular_SeqSBAIJ, 1399 MatRestoreRowUpperTriangular_SeqSBAIJ, 1400 /*109*/0, 1401 0, 1402 0, 1403 0, 1404 MatMissingDiagonal_SeqSBAIJ, 1405 /*114*/0, 1406 0, 1407 0, 1408 0, 1409 0, 1410 /*119*/0, 1411 0, 1412 0 1413 }; 1414 1415 EXTERN_C_BEGIN 1416 #undef __FUNCT__ 1417 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1418 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1419 { 1420 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1421 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1422 PetscErrorCode ierr; 1423 1424 PetscFunctionBegin; 1425 if (aij->nonew != 1) { 1426 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1427 } 1428 1429 /* allocate space for values if not already there */ 1430 if (!aij->saved_values) { 1431 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1432 } 1433 1434 /* copy values over */ 1435 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 EXTERN_C_END 1439 1440 EXTERN_C_BEGIN 1441 #undef __FUNCT__ 1442 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1443 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1444 { 1445 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1446 PetscErrorCode ierr; 1447 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1448 1449 PetscFunctionBegin; 1450 if (aij->nonew != 1) { 1451 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1452 } 1453 if (!aij->saved_values) { 1454 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1455 } 1456 1457 /* copy values over */ 1458 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1459 PetscFunctionReturn(0); 1460 } 1461 EXTERN_C_END 1462 1463 EXTERN_C_BEGIN 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1466 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1467 { 1468 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1469 PetscErrorCode ierr; 1470 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1471 PetscTruth skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1472 1473 PetscFunctionBegin; 1474 B->preallocated = PETSC_TRUE; 1475 if (bs < 0) { 1476 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1477 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1478 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1479 bs = PetscAbs(bs); 1480 } 1481 if (nnz && newbs != bs) { 1482 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1483 } 1484 bs = newbs; 1485 1486 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1487 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1488 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1489 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1490 1491 mbs = B->rmap->N/bs; 1492 bs2 = bs*bs; 1493 1494 if (mbs*bs != B->rmap->N) { 1495 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1496 } 1497 1498 if (nz == MAT_SKIP_ALLOCATION) { 1499 skipallocation = PETSC_TRUE; 1500 nz = 0; 1501 } 1502 1503 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1504 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1505 if (nnz) { 1506 for (i=0; i<mbs; i++) { 1507 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1508 if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],mbs); 1509 } 1510 } 1511 1512 B->ops->mult = MatMult_SeqSBAIJ_N; 1513 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1514 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1515 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1516 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1517 if (!flg) { 1518 switch (bs) { 1519 case 1: 1520 B->ops->mult = MatMult_SeqSBAIJ_1; 1521 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1522 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1523 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1524 break; 1525 case 2: 1526 B->ops->mult = MatMult_SeqSBAIJ_2; 1527 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1528 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1529 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1530 break; 1531 case 3: 1532 B->ops->mult = MatMult_SeqSBAIJ_3; 1533 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1534 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1535 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1536 break; 1537 case 4: 1538 B->ops->mult = MatMult_SeqSBAIJ_4; 1539 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1540 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1541 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1542 break; 1543 case 5: 1544 B->ops->mult = MatMult_SeqSBAIJ_5; 1545 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1546 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1547 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1548 break; 1549 case 6: 1550 B->ops->mult = MatMult_SeqSBAIJ_6; 1551 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1552 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1553 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1554 break; 1555 case 7: 1556 B->ops->mult = MatMult_SeqSBAIJ_7; 1557 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1558 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1559 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1560 break; 1561 } 1562 } 1563 1564 b->mbs = mbs; 1565 b->nbs = mbs; 1566 if (!skipallocation) { 1567 if (!b->imax) { 1568 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1569 b->free_imax_ilen = PETSC_TRUE; 1570 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1571 } 1572 if (!nnz) { 1573 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1574 else if (nz <= 0) nz = 1; 1575 for (i=0; i<mbs; i++) { 1576 b->imax[i] = nz; 1577 } 1578 nz = nz*mbs; /* total nz */ 1579 } else { 1580 nz = 0; 1581 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1582 } 1583 /* b->ilen will count nonzeros in each block row so far. */ 1584 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1585 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1586 1587 /* allocate the matrix space */ 1588 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1589 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1590 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1591 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1592 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1593 b->singlemalloc = PETSC_TRUE; 1594 1595 /* pointer to beginning of each row */ 1596 b->i[0] = 0; 1597 for (i=1; i<mbs+1; i++) { 1598 b->i[i] = b->i[i-1] + b->imax[i-1]; 1599 } 1600 b->free_a = PETSC_TRUE; 1601 b->free_ij = PETSC_TRUE; 1602 } else { 1603 b->free_a = PETSC_FALSE; 1604 b->free_ij = PETSC_FALSE; 1605 } 1606 1607 B->rmap->bs = bs; 1608 b->bs2 = bs2; 1609 b->nz = 0; 1610 b->maxnz = nz*bs2; 1611 1612 b->inew = 0; 1613 b->jnew = 0; 1614 b->anew = 0; 1615 b->a2anew = 0; 1616 b->permute = PETSC_FALSE; 1617 PetscFunctionReturn(0); 1618 } 1619 EXTERN_C_END 1620 1621 /* 1622 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1623 */ 1624 #undef __FUNCT__ 1625 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1626 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscTruth natural) 1627 { 1628 PetscErrorCode ierr; 1629 PetscTruth flg = PETSC_FALSE; 1630 PetscInt bs = B->rmap->bs; 1631 1632 PetscFunctionBegin; 1633 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1634 if (flg) bs = 8; 1635 1636 if (!natural) { 1637 switch (bs) { 1638 case 1: 1639 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1640 break; 1641 case 2: 1642 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1643 break; 1644 case 3: 1645 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1646 break; 1647 case 4: 1648 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1649 break; 1650 case 5: 1651 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1652 break; 1653 case 6: 1654 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1655 break; 1656 case 7: 1657 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1658 break; 1659 default: 1660 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1661 break; 1662 } 1663 } else { 1664 switch (bs) { 1665 case 1: 1666 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1667 break; 1668 case 2: 1669 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1670 break; 1671 case 3: 1672 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1673 break; 1674 case 4: 1675 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1676 break; 1677 case 5: 1678 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1679 break; 1680 case 6: 1681 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1682 break; 1683 case 7: 1684 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1685 break; 1686 default: 1687 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1688 break; 1689 } 1690 } 1691 PetscFunctionReturn(0); 1692 } 1693 1694 EXTERN_C_BEGIN 1695 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1696 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1697 EXTERN_C_END 1698 1699 1700 EXTERN_C_BEGIN 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1703 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1704 { 1705 PetscInt n = A->rmap->n; 1706 PetscErrorCode ierr; 1707 1708 PetscFunctionBegin; 1709 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1710 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1711 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1712 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1713 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1714 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1715 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1716 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1717 (*B)->factor = ftype; 1718 PetscFunctionReturn(0); 1719 } 1720 EXTERN_C_END 1721 1722 EXTERN_C_BEGIN 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1725 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1726 { 1727 PetscFunctionBegin; 1728 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1729 *flg = PETSC_TRUE; 1730 } else { 1731 *flg = PETSC_FALSE; 1732 } 1733 PetscFunctionReturn(0); 1734 } 1735 EXTERN_C_END 1736 1737 EXTERN_C_BEGIN 1738 #if defined(PETSC_HAVE_MUMPS) 1739 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1740 #endif 1741 #if defined(PETSC_HAVE_SPOOLES) 1742 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1743 #endif 1744 #if defined(PETSC_HAVE_PASTIX) 1745 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1746 #endif 1747 EXTERN_C_END 1748 1749 /*MC 1750 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1751 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1752 1753 Options Database Keys: 1754 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1755 1756 Level: beginner 1757 1758 .seealso: MatCreateSeqSBAIJ 1759 M*/ 1760 1761 EXTERN_C_BEGIN 1762 #undef __FUNCT__ 1763 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1764 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1765 { 1766 Mat_SeqSBAIJ *b; 1767 PetscErrorCode ierr; 1768 PetscMPIInt size; 1769 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1770 1771 PetscFunctionBegin; 1772 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1773 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1774 1775 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1776 B->data = (void*)b; 1777 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1778 B->ops->destroy = MatDestroy_SeqSBAIJ; 1779 B->ops->view = MatView_SeqSBAIJ; 1780 B->mapping = 0; 1781 b->row = 0; 1782 b->icol = 0; 1783 b->reallocs = 0; 1784 b->saved_values = 0; 1785 b->inode.limit = 5; 1786 b->inode.max_limit = 5; 1787 1788 b->roworiented = PETSC_TRUE; 1789 b->nonew = 0; 1790 b->diag = 0; 1791 b->solve_work = 0; 1792 b->mult_work = 0; 1793 B->spptr = 0; 1794 b->keepnonzeropattern = PETSC_FALSE; 1795 b->xtoy = 0; 1796 b->XtoY = 0; 1797 1798 b->inew = 0; 1799 b->jnew = 0; 1800 b->anew = 0; 1801 b->a2anew = 0; 1802 b->permute = PETSC_FALSE; 1803 1804 b->ignore_ltriangular = PETSC_FALSE; 1805 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1806 1807 b->getrow_utriangular = PETSC_FALSE; 1808 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1809 1810 #if defined(PETSC_HAVE_PASTIX) 1811 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1812 "MatGetFactor_seqsbaij_pastix", 1813 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1814 #endif 1815 #if defined(PETSC_HAVE_SPOOLES) 1816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1817 "MatGetFactor_seqsbaij_spooles", 1818 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1819 #endif 1820 #if defined(PETSC_HAVE_MUMPS) 1821 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1822 "MatGetFactor_seqsbaij_mumps", 1823 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1824 #endif 1825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1826 "MatGetFactorAvailable_seqsbaij_petsc", 1827 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1828 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1829 "MatGetFactor_seqsbaij_petsc", 1830 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1831 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1832 "MatStoreValues_SeqSBAIJ", 1833 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1834 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1835 "MatRetrieveValues_SeqSBAIJ", 1836 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1837 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1838 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1839 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1840 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1841 "MatConvert_SeqSBAIJ_SeqAIJ", 1842 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1843 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1844 "MatConvert_SeqSBAIJ_SeqBAIJ", 1845 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1846 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1847 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1848 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1849 1850 B->symmetric = PETSC_TRUE; 1851 B->structurally_symmetric = PETSC_TRUE; 1852 B->symmetric_set = PETSC_TRUE; 1853 B->structurally_symmetric_set = PETSC_TRUE; 1854 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1855 1856 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1857 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1858 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1859 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1860 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1861 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",PETSC_NULL,b->inode.limit,&b->inode.limit,PETSC_NULL);CHKERRQ(ierr); 1862 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1863 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1864 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1865 1866 PetscFunctionReturn(0); 1867 } 1868 EXTERN_C_END 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1872 /*@C 1873 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1874 compressed row) format. For good matrix assembly performance the 1875 user should preallocate the matrix storage by setting the parameter nz 1876 (or the array nnz). By setting these parameters accurately, performance 1877 during matrix assembly can be increased by more than a factor of 50. 1878 1879 Collective on Mat 1880 1881 Input Parameters: 1882 + A - the symmetric matrix 1883 . bs - size of block 1884 . nz - number of block nonzeros per block row (same for all rows) 1885 - nnz - array containing the number of block nonzeros in the upper triangular plus 1886 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1887 1888 Options Database Keys: 1889 . -mat_no_unroll - uses code that does not unroll the loops in the 1890 block calculations (much slower) 1891 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1892 1893 Level: intermediate 1894 1895 Notes: 1896 Specify the preallocated storage with either nz or nnz (not both). 1897 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1898 allocation. For additional details, see the users manual chapter on 1899 matrices. 1900 1901 You can call MatGetInfo() to get information on how effective the preallocation was; 1902 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1903 You can also run with the option -info and look for messages with the string 1904 malloc in them to see if additional memory allocation was needed. 1905 1906 If the nnz parameter is given then the nz parameter is ignored 1907 1908 1909 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1910 @*/ 1911 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1912 { 1913 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1914 1915 PetscFunctionBegin; 1916 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1917 if (f) { 1918 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1919 } 1920 PetscFunctionReturn(0); 1921 } 1922 1923 #undef __FUNCT__ 1924 #define __FUNCT__ "MatCreateSeqSBAIJ" 1925 /*@C 1926 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1927 compressed row) format. For good matrix assembly performance the 1928 user should preallocate the matrix storage by setting the parameter nz 1929 (or the array nnz). By setting these parameters accurately, performance 1930 during matrix assembly can be increased by more than a factor of 50. 1931 1932 Collective on MPI_Comm 1933 1934 Input Parameters: 1935 + comm - MPI communicator, set to PETSC_COMM_SELF 1936 . bs - size of block 1937 . m - number of rows, or number of columns 1938 . nz - number of block nonzeros per block row (same for all rows) 1939 - nnz - array containing the number of block nonzeros in the upper triangular plus 1940 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1941 1942 Output Parameter: 1943 . A - the symmetric matrix 1944 1945 Options Database Keys: 1946 . -mat_no_unroll - uses code that does not unroll the loops in the 1947 block calculations (much slower) 1948 . -mat_block_size - size of the blocks to use 1949 1950 Level: intermediate 1951 1952 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1953 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1954 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1955 1956 Notes: 1957 The number of rows and columns must be divisible by blocksize. 1958 This matrix type does not support complex Hermitian operation. 1959 1960 Specify the preallocated storage with either nz or nnz (not both). 1961 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1962 allocation. For additional details, see the users manual chapter on 1963 matrices. 1964 1965 If the nnz parameter is given then the nz parameter is ignored 1966 1967 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1968 @*/ 1969 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1970 { 1971 PetscErrorCode ierr; 1972 1973 PetscFunctionBegin; 1974 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1975 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1976 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1977 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1978 PetscFunctionReturn(0); 1979 } 1980 1981 #undef __FUNCT__ 1982 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1983 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1984 { 1985 Mat C; 1986 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1987 PetscErrorCode ierr; 1988 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1989 1990 PetscFunctionBegin; 1991 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1992 1993 *B = 0; 1994 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1995 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1996 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1997 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1998 c = (Mat_SeqSBAIJ*)C->data; 1999 2000 C->preallocated = PETSC_TRUE; 2001 C->factor = A->factor; 2002 c->row = 0; 2003 c->icol = 0; 2004 c->saved_values = 0; 2005 c->keepnonzeropattern = a->keepnonzeropattern; 2006 C->assembled = PETSC_TRUE; 2007 2008 ierr = PetscLayoutCopy(A->rmap,&C->rmap);CHKERRQ(ierr); 2009 ierr = PetscLayoutCopy(A->cmap,&C->cmap);CHKERRQ(ierr); 2010 c->bs2 = a->bs2; 2011 c->mbs = a->mbs; 2012 c->nbs = a->nbs; 2013 2014 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2015 c->imax = a->imax; 2016 c->ilen = a->ilen; 2017 c->free_imax_ilen = PETSC_FALSE; 2018 } else { 2019 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2020 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2021 for (i=0; i<mbs; i++) { 2022 c->imax[i] = a->imax[i]; 2023 c->ilen[i] = a->ilen[i]; 2024 } 2025 c->free_imax_ilen = PETSC_TRUE; 2026 } 2027 2028 /* allocate the matrix space */ 2029 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2030 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2031 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2032 c->singlemalloc = PETSC_FALSE; 2033 c->free_ij = PETSC_FALSE; 2034 c->parent = A; 2035 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2036 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2037 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2038 } else { 2039 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2040 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2041 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2042 c->singlemalloc = PETSC_TRUE; 2043 c->free_ij = PETSC_TRUE; 2044 } 2045 if (mbs > 0) { 2046 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2047 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2048 } 2049 if (cpvalues == MAT_COPY_VALUES) { 2050 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2051 } else { 2052 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2053 } 2054 if (a->jshort) { 2055 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2056 c->jshort = a->jshort; 2057 c->free_jshort = PETSC_FALSE; 2058 } else { 2059 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2060 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2061 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2062 c->free_jshort = PETSC_TRUE; 2063 } 2064 } 2065 } 2066 2067 c->roworiented = a->roworiented; 2068 c->nonew = a->nonew; 2069 2070 if (a->diag) { 2071 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2072 c->diag = a->diag; 2073 c->free_diag = PETSC_FALSE; 2074 } else { 2075 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2076 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2077 for (i=0; i<mbs; i++) { 2078 c->diag[i] = a->diag[i]; 2079 } 2080 c->free_diag = PETSC_TRUE; 2081 } 2082 } else c->diag = 0; 2083 c->nz = a->nz; 2084 c->maxnz = a->maxnz; 2085 c->solve_work = 0; 2086 c->mult_work = 0; 2087 c->free_a = PETSC_TRUE; 2088 *B = C; 2089 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 #undef __FUNCT__ 2094 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2095 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 2096 { 2097 Mat_SeqSBAIJ *a; 2098 Mat B; 2099 PetscErrorCode ierr; 2100 int fd; 2101 PetscMPIInt size; 2102 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2103 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2104 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2105 PetscInt *masked,nmask,tmp,bs2,ishift; 2106 PetscScalar *aa; 2107 MPI_Comm comm = ((PetscObject)viewer)->comm; 2108 2109 PetscFunctionBegin; 2110 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2111 bs2 = bs*bs; 2112 2113 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2114 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 2115 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2116 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2117 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2118 M = header[1]; N = header[2]; nz = header[3]; 2119 2120 if (header[3] < 0) { 2121 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2122 } 2123 2124 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 2125 2126 /* 2127 This code adds extra rows to make sure the number of rows is 2128 divisible by the blocksize 2129 */ 2130 mbs = M/bs; 2131 extra_rows = bs - M + bs*(mbs); 2132 if (extra_rows == bs) extra_rows = 0; 2133 else mbs++; 2134 if (extra_rows) { 2135 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2136 } 2137 2138 /* read in row lengths */ 2139 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2140 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2141 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2142 2143 /* read in column indices */ 2144 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2145 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2146 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2147 2148 /* loop over row lengths determining block row lengths */ 2149 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2150 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2151 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2152 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2153 rowcount = 0; 2154 nzcount = 0; 2155 for (i=0; i<mbs; i++) { 2156 nmask = 0; 2157 for (j=0; j<bs; j++) { 2158 kmax = rowlengths[rowcount]; 2159 for (k=0; k<kmax; k++) { 2160 tmp = jj[nzcount++]/bs; /* block col. index */ 2161 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2162 } 2163 rowcount++; 2164 } 2165 s_browlengths[i] += nmask; 2166 2167 /* zero out the mask elements we set */ 2168 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2169 } 2170 2171 /* create our matrix */ 2172 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2173 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2174 ierr = MatSetType(B,type);CHKERRQ(ierr); 2175 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2176 a = (Mat_SeqSBAIJ*)B->data; 2177 2178 /* set matrix "i" values */ 2179 a->i[0] = 0; 2180 for (i=1; i<= mbs; i++) { 2181 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2182 a->ilen[i-1] = s_browlengths[i-1]; 2183 } 2184 a->nz = a->i[mbs]; 2185 2186 /* read in nonzero values */ 2187 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2188 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2189 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2190 2191 /* set "a" and "j" values into matrix */ 2192 nzcount = 0; jcount = 0; 2193 for (i=0; i<mbs; i++) { 2194 nzcountb = nzcount; 2195 nmask = 0; 2196 for (j=0; j<bs; j++) { 2197 kmax = rowlengths[i*bs+j]; 2198 for (k=0; k<kmax; k++) { 2199 tmp = jj[nzcount++]/bs; /* block col. index */ 2200 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2201 } 2202 } 2203 /* sort the masked values */ 2204 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2205 2206 /* set "j" values into matrix */ 2207 maskcount = 1; 2208 for (j=0; j<nmask; j++) { 2209 a->j[jcount++] = masked[j]; 2210 mask[masked[j]] = maskcount++; 2211 } 2212 2213 /* set "a" values into matrix */ 2214 ishift = bs2*a->i[i]; 2215 for (j=0; j<bs; j++) { 2216 kmax = rowlengths[i*bs+j]; 2217 for (k=0; k<kmax; k++) { 2218 tmp = jj[nzcountb]/bs ; /* block col. index */ 2219 if (tmp >= i){ 2220 block = mask[tmp] - 1; 2221 point = jj[nzcountb] - bs*tmp; 2222 idx = ishift + bs2*block + j + bs*point; 2223 a->a[idx] = aa[nzcountb]; 2224 } 2225 nzcountb++; 2226 } 2227 } 2228 /* zero out the mask elements we set */ 2229 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2230 } 2231 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2232 2233 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2234 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2235 ierr = PetscFree(aa);CHKERRQ(ierr); 2236 ierr = PetscFree(jj);CHKERRQ(ierr); 2237 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2238 2239 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2240 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2241 ierr = MatView_Private(B);CHKERRQ(ierr); 2242 *A = B; 2243 PetscFunctionReturn(0); 2244 } 2245 2246 #undef __FUNCT__ 2247 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2248 /*@ 2249 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2250 (upper triangular entries in CSR format) provided by the user. 2251 2252 Collective on MPI_Comm 2253 2254 Input Parameters: 2255 + comm - must be an MPI communicator of size 1 2256 . bs - size of block 2257 . m - number of rows 2258 . n - number of columns 2259 . i - row indices 2260 . j - column indices 2261 - a - matrix values 2262 2263 Output Parameter: 2264 . mat - the matrix 2265 2266 Level: intermediate 2267 2268 Notes: 2269 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2270 once the matrix is destroyed 2271 2272 You cannot set new nonzero locations into this matrix, that will generate an error. 2273 2274 The i and j indices are 0 based 2275 2276 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2277 2278 @*/ 2279 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2280 { 2281 PetscErrorCode ierr; 2282 PetscInt ii; 2283 Mat_SeqSBAIJ *sbaij; 2284 2285 PetscFunctionBegin; 2286 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2287 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2288 2289 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2290 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2291 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2292 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2293 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2294 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2295 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2296 2297 sbaij->i = i; 2298 sbaij->j = j; 2299 sbaij->a = a; 2300 sbaij->singlemalloc = PETSC_FALSE; 2301 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2302 sbaij->free_a = PETSC_FALSE; 2303 sbaij->free_ij = PETSC_FALSE; 2304 2305 for (ii=0; ii<m; ii++) { 2306 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2307 #if defined(PETSC_USE_DEBUG) 2308 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2309 #endif 2310 } 2311 #if defined(PETSC_USE_DEBUG) 2312 for (ii=0; ii<sbaij->i[m]; ii++) { 2313 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2314 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2315 } 2316 #endif 2317 2318 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2319 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 2324 2325 2326 2327