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(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 defined(USESHORT) 194 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort; 195 #else 196 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian; 197 #endif 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 533 PetscFunctionBegin; 534 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 535 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 536 if (iascii){ 537 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 538 } else if (isdraw) { 539 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 540 } else { 541 Mat B; 542 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 543 ierr = MatView(B,viewer);CHKERRQ(ierr); 544 ierr = MatDestroy(B);CHKERRQ(ierr); 545 } 546 PetscFunctionReturn(0); 547 } 548 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 552 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 553 { 554 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 555 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 556 PetscInt *ai = a->i,*ailen = a->ilen; 557 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 558 MatScalar *ap,*aa = a->a; 559 560 PetscFunctionBegin; 561 for (k=0; k<m; k++) { /* loop over rows */ 562 row = im[k]; brow = row/bs; 563 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 564 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 565 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 566 nrow = ailen[brow]; 567 for (l=0; l<n; l++) { /* loop over columns */ 568 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 569 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 570 col = in[l] ; 571 bcol = col/bs; 572 cidx = col%bs; 573 ridx = row%bs; 574 high = nrow; 575 low = 0; /* assume unsorted */ 576 while (high-low > 5) { 577 t = (low+high)/2; 578 if (rp[t] > bcol) high = t; 579 else low = t; 580 } 581 for (i=low; i<high; i++) { 582 if (rp[i] > bcol) break; 583 if (rp[i] == bcol) { 584 *v++ = ap[bs2*i+bs*cidx+ridx]; 585 goto finished; 586 } 587 } 588 *v++ = 0.0; 589 finished:; 590 } 591 } 592 PetscFunctionReturn(0); 593 } 594 595 596 #undef __FUNCT__ 597 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 598 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 599 { 600 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 601 PetscErrorCode ierr; 602 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 603 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 604 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 605 PetscTruth roworiented=a->roworiented; 606 const PetscScalar *value = v; 607 MatScalar *ap,*aa = a->a,*bap; 608 609 PetscFunctionBegin; 610 if (roworiented) { 611 stepval = (n-1)*bs; 612 } else { 613 stepval = (m-1)*bs; 614 } 615 for (k=0; k<m; k++) { /* loop over added rows */ 616 row = im[k]; 617 if (row < 0) continue; 618 #if defined(PETSC_USE_DEBUG) 619 if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 620 #endif 621 rp = aj + ai[row]; 622 ap = aa + bs2*ai[row]; 623 rmax = imax[row]; 624 nrow = ailen[row]; 625 low = 0; 626 high = nrow; 627 for (l=0; l<n; l++) { /* loop over added columns */ 628 if (in[l] < 0) continue; 629 col = in[l]; 630 #if defined(PETSC_USE_DEBUG) 631 if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 632 #endif 633 if (col < row) continue; /* ignore lower triangular block */ 634 if (roworiented) { 635 value = v + k*(stepval+bs)*bs + l*bs; 636 } else { 637 value = v + l*(stepval+bs)*bs + k*bs; 638 } 639 if (col <= lastcol) low = 0; else high = nrow; 640 lastcol = col; 641 while (high-low > 7) { 642 t = (low+high)/2; 643 if (rp[t] > col) high = t; 644 else low = t; 645 } 646 for (i=low; i<high; i++) { 647 if (rp[i] > col) break; 648 if (rp[i] == col) { 649 bap = ap + bs2*i; 650 if (roworiented) { 651 if (is == ADD_VALUES) { 652 for (ii=0; ii<bs; ii++,value+=stepval) { 653 for (jj=ii; jj<bs2; jj+=bs) { 654 bap[jj] += *value++; 655 } 656 } 657 } else { 658 for (ii=0; ii<bs; ii++,value+=stepval) { 659 for (jj=ii; jj<bs2; jj+=bs) { 660 bap[jj] = *value++; 661 } 662 } 663 } 664 } else { 665 if (is == ADD_VALUES) { 666 for (ii=0; ii<bs; ii++,value+=stepval) { 667 for (jj=0; jj<bs; jj++) { 668 *bap++ += *value++; 669 } 670 } 671 } else { 672 for (ii=0; ii<bs; ii++,value+=stepval) { 673 for (jj=0; jj<bs; jj++) { 674 *bap++ = *value++; 675 } 676 } 677 } 678 } 679 goto noinsert2; 680 } 681 } 682 if (nonew == 1) goto noinsert2; 683 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 684 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 685 N = nrow++ - 1; high++; 686 /* shift up all the later entries in this row */ 687 for (ii=N; ii>=i; ii--) { 688 rp[ii+1] = rp[ii]; 689 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 690 } 691 if (N >= i) { 692 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 693 } 694 rp[i] = col; 695 bap = ap + bs2*i; 696 if (roworiented) { 697 for (ii=0; ii<bs; ii++,value+=stepval) { 698 for (jj=ii; jj<bs2; jj+=bs) { 699 bap[jj] = *value++; 700 } 701 } 702 } else { 703 for (ii=0; ii<bs; ii++,value+=stepval) { 704 for (jj=0; jj<bs; jj++) { 705 *bap++ = *value++; 706 } 707 } 708 } 709 noinsert2:; 710 low = i; 711 } 712 ailen[row] = nrow; 713 } 714 PetscFunctionReturn(0); 715 } 716 717 /* 718 This is not yet used 719 */ 720 #undef __FUNCT__ 721 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 722 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 723 { 724 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 725 PetscErrorCode ierr; 726 const PetscInt *ai = a->i, *aj = a->j,*cols; 727 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 728 PetscTruth flag; 729 730 PetscFunctionBegin; 731 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 732 while (i < m){ 733 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 734 /* Limits the number of elements in a node to 'a->inode.limit' */ 735 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 736 nzy = ai[j+1] - ai[j]; 737 if (nzy != (nzx - j + i)) break; 738 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 739 if (!flag) break; 740 } 741 ns[node_count++] = blk_size; 742 i = j; 743 } 744 if (!a->inode.size && m && node_count > .9*m) { 745 ierr = PetscFree(ns);CHKERRQ(ierr); 746 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 747 } else { 748 a->inode.node_count = node_count; 749 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 750 ierr = PetscLogObjectMemory(A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 751 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 752 ierr = PetscFree(ns);CHKERRQ(ierr); 753 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 754 755 /* count collections of adjacent columns in each inode */ 756 row = 0; 757 cnt = 0; 758 for (i=0; i<node_count; i++) { 759 cols = aj + ai[row] + a->inode.size[i]; 760 nz = ai[row+1] - ai[row] - a->inode.size[i]; 761 for (j=1; j<nz; j++) { 762 if (cols[j] != cols[j-1]+1) { 763 cnt++; 764 } 765 } 766 cnt++; 767 row += a->inode.size[i]; 768 } 769 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 770 cnt = 0; 771 row = 0; 772 for (i=0; i<node_count; i++) { 773 cols = aj + ai[row] + a->inode.size[i]; 774 CHKMEMQ; 775 counts[2*cnt] = cols[0]; 776 CHKMEMQ; 777 nz = ai[row+1] - ai[row] - a->inode.size[i]; 778 cnt2 = 1; 779 for (j=1; j<nz; j++) { 780 if (cols[j] != cols[j-1]+1) { 781 CHKMEMQ; 782 counts[2*(cnt++)+1] = cnt2; 783 counts[2*cnt] = cols[j]; 784 CHKMEMQ; 785 cnt2 = 1; 786 } else cnt2++; 787 } 788 CHKMEMQ; 789 counts[2*(cnt++)+1] = cnt2; 790 CHKMEMQ; 791 row += a->inode.size[i]; 792 } 793 ierr = PetscIntView(2*cnt,counts,0); 794 } 795 PetscFunctionReturn(0); 796 } 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 800 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 801 { 802 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 803 PetscErrorCode ierr; 804 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 805 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 806 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 807 MatScalar *aa = a->a,*ap; 808 809 PetscFunctionBegin; 810 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 811 812 if (m) rmax = ailen[0]; 813 for (i=1; i<mbs; i++) { 814 /* move each row back by the amount of empty slots (fshift) before it*/ 815 fshift += imax[i-1] - ailen[i-1]; 816 rmax = PetscMax(rmax,ailen[i]); 817 if (fshift) { 818 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 819 N = ailen[i]; 820 for (j=0; j<N; j++) { 821 ip[j-fshift] = ip[j]; 822 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 823 } 824 } 825 ai[i] = ai[i-1] + ailen[i-1]; 826 } 827 if (mbs) { 828 fshift += imax[mbs-1] - ailen[mbs-1]; 829 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 830 } 831 /* reset ilen and imax for each row */ 832 for (i=0; i<mbs; i++) { 833 ailen[i] = imax[i] = ai[i+1] - ai[i]; 834 } 835 a->nz = ai[mbs]; 836 837 /* diagonals may have moved, reset it */ 838 if (a->diag) { 839 ierr = PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 840 } 841 if (fshift && a->nounused == -1) { 842 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); 843 } 844 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); 845 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 846 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 847 a->reallocs = 0; 848 A->info.nz_unneeded = (PetscReal)fshift*bs2; 849 a->idiagvalid = PETSC_FALSE; 850 851 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 852 if (!a->jshort) { 853 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 854 ierr = PetscLogObjectMemory(A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 855 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 856 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 857 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 858 a->free_jshort = PETSC_TRUE; 859 } 860 } 861 PetscFunctionReturn(0); 862 } 863 864 /* 865 This function returns an array of flags which indicate the locations of contiguous 866 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 867 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 868 Assume: sizes should be long enough to hold all the values. 869 */ 870 #undef __FUNCT__ 871 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 872 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 873 { 874 PetscInt i,j,k,row; 875 PetscTruth flg; 876 877 PetscFunctionBegin; 878 for (i=0,j=0; i<n; j++) { 879 row = idx[i]; 880 if (row%bs!=0) { /* Not the begining of a block */ 881 sizes[j] = 1; 882 i++; 883 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 884 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 885 i++; 886 } else { /* Begining of the block, so check if the complete block exists */ 887 flg = PETSC_TRUE; 888 for (k=1; k<bs; k++) { 889 if (row+k != idx[i+k]) { /* break in the block */ 890 flg = PETSC_FALSE; 891 break; 892 } 893 } 894 if (flg) { /* No break in the bs */ 895 sizes[j] = bs; 896 i+= bs; 897 } else { 898 sizes[j] = 1; 899 i++; 900 } 901 } 902 } 903 *bs_max = j; 904 PetscFunctionReturn(0); 905 } 906 907 908 /* Only add/insert a(i,j) with i<=j (blocks). 909 Any a(i,j) with i>j input by user is ingored. 910 */ 911 912 #undef __FUNCT__ 913 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 914 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 915 { 916 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 917 PetscErrorCode ierr; 918 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 919 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 920 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 921 PetscInt ridx,cidx,bs2=a->bs2; 922 MatScalar *ap,value,*aa=a->a,*bap; 923 924 PetscFunctionBegin; 925 for (k=0; k<m; k++) { /* loop over added rows */ 926 row = im[k]; /* row number */ 927 brow = row/bs; /* block row number */ 928 if (row < 0) continue; 929 #if defined(PETSC_USE_DEBUG) 930 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 931 #endif 932 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 933 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 934 rmax = imax[brow]; /* maximum space allocated for this row */ 935 nrow = ailen[brow]; /* actual length of this row */ 936 low = 0; 937 938 for (l=0; l<n; l++) { /* loop over added columns */ 939 if (in[l] < 0) continue; 940 #if defined(PETSC_USE_DEBUG) 941 if (in[l] >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->N-1); 942 #endif 943 col = in[l]; 944 bcol = col/bs; /* block col number */ 945 946 if (brow > bcol) { 947 if (a->ignore_ltriangular){ 948 continue; /* ignore lower triangular values */ 949 } else { 950 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)"); 951 } 952 } 953 954 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 955 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 956 /* element value a(k,l) */ 957 if (roworiented) { 958 value = v[l + k*n]; 959 } else { 960 value = v[k + l*m]; 961 } 962 963 /* move pointer bap to a(k,l) quickly and add/insert value */ 964 if (col <= lastcol) low = 0; high = nrow; 965 lastcol = col; 966 while (high-low > 7) { 967 t = (low+high)/2; 968 if (rp[t] > bcol) high = t; 969 else low = t; 970 } 971 for (i=low; i<high; i++) { 972 if (rp[i] > bcol) break; 973 if (rp[i] == bcol) { 974 bap = ap + bs2*i + bs*cidx + ridx; 975 if (is == ADD_VALUES) *bap += value; 976 else *bap = value; 977 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 978 if (brow == bcol && ridx < cidx){ 979 bap = ap + bs2*i + bs*ridx + cidx; 980 if (is == ADD_VALUES) *bap += value; 981 else *bap = value; 982 } 983 goto noinsert1; 984 } 985 } 986 987 if (nonew == 1) goto noinsert1; 988 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 989 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 990 991 N = nrow++ - 1; high++; 992 /* shift up all the later entries in this row */ 993 for (ii=N; ii>=i; ii--) { 994 rp[ii+1] = rp[ii]; 995 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 996 } 997 if (N>=i) { 998 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 999 } 1000 rp[i] = bcol; 1001 ap[bs2*i + bs*cidx + ridx] = value; 1002 noinsert1:; 1003 low = i; 1004 } 1005 } /* end of loop over added columns */ 1006 ailen[brow] = nrow; 1007 } /* end of loop over added rows */ 1008 PetscFunctionReturn(0); 1009 } 1010 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 1013 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1014 { 1015 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1016 Mat outA; 1017 PetscErrorCode ierr; 1018 PetscTruth row_identity; 1019 1020 PetscFunctionBegin; 1021 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1022 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1023 if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported"); 1024 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()! */ 1025 1026 outA = inA; 1027 inA->factor = MAT_FACTOR_ICC; 1028 1029 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1030 ierr = MatSeqSBAIJSetNumericFactorization(inA,row_identity);CHKERRQ(ierr); 1031 1032 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1033 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr); } 1034 a->row = row; 1035 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1036 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr); } 1037 a->col = row; 1038 1039 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1040 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1041 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1042 1043 if (!a->solve_work) { 1044 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1045 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1046 } 1047 1048 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1049 PetscFunctionReturn(0); 1050 } 1051 1052 EXTERN_C_BEGIN 1053 #undef __FUNCT__ 1054 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1055 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1056 { 1057 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1058 PetscInt i,nz,n; 1059 1060 PetscFunctionBegin; 1061 nz = baij->maxnz; 1062 n = mat->cmap->n; 1063 for (i=0; i<nz; i++) { 1064 baij->j[i] = indices[i]; 1065 } 1066 baij->nz = nz; 1067 for (i=0; i<n; i++) { 1068 baij->ilen[i] = baij->imax[i]; 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 EXTERN_C_END 1073 1074 #undef __FUNCT__ 1075 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1076 /*@ 1077 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1078 in the matrix. 1079 1080 Input Parameters: 1081 + mat - the SeqSBAIJ matrix 1082 - indices - the column indices 1083 1084 Level: advanced 1085 1086 Notes: 1087 This can be called if you have precomputed the nonzero structure of the 1088 matrix and want to provide it to the matrix object to improve the performance 1089 of the MatSetValues() operation. 1090 1091 You MUST have set the correct numbers of nonzeros per row in the call to 1092 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1093 1094 MUST be called before any calls to MatSetValues() 1095 1096 .seealso: MatCreateSeqSBAIJ 1097 @*/ 1098 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1099 { 1100 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1101 1102 PetscFunctionBegin; 1103 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1104 PetscValidPointer(indices,2); 1105 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1106 if (f) { 1107 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1108 } else { 1109 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1110 } 1111 PetscFunctionReturn(0); 1112 } 1113 1114 #undef __FUNCT__ 1115 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1116 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1117 { 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 /* If the two matrices have the same copy implementation, use fast copy. */ 1122 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1123 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1124 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1125 1126 if (a->i[A->rmap->N] != b->i[B->rmap->N]) { 1127 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1128 } 1129 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1130 } else { 1131 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1132 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1133 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1134 } 1135 PetscFunctionReturn(0); 1136 } 1137 1138 #undef __FUNCT__ 1139 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1140 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1141 { 1142 PetscErrorCode ierr; 1143 1144 PetscFunctionBegin; 1145 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1151 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1152 { 1153 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1154 PetscFunctionBegin; 1155 *array = a->a; 1156 PetscFunctionReturn(0); 1157 } 1158 1159 #undef __FUNCT__ 1160 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1161 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1162 { 1163 PetscFunctionBegin; 1164 PetscFunctionReturn(0); 1165 } 1166 1167 #include "petscblaslapack.h" 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1170 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1171 { 1172 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1173 PetscErrorCode ierr; 1174 PetscInt i,bs=Y->rmap->bs,bs2,j; 1175 PetscBLASInt one = 1,bnz = PetscBLASIntCast(x->nz); 1176 1177 PetscFunctionBegin; 1178 if (str == SAME_NONZERO_PATTERN) { 1179 PetscScalar alpha = a; 1180 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1181 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1182 if (y->xtoy && y->XtoY != X) { 1183 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1184 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1185 } 1186 if (!y->xtoy) { /* get xtoy */ 1187 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1188 y->XtoY = X; 1189 } 1190 bs2 = bs*bs; 1191 for (i=0; i<x->nz; i++) { 1192 j = 0; 1193 while (j < bs2){ 1194 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1195 j++; 1196 } 1197 } 1198 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); 1199 } else { 1200 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1201 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1202 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1203 } 1204 PetscFunctionReturn(0); 1205 } 1206 1207 #undef __FUNCT__ 1208 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1209 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1210 { 1211 PetscFunctionBegin; 1212 *flg = PETSC_TRUE; 1213 PetscFunctionReturn(0); 1214 } 1215 1216 #undef __FUNCT__ 1217 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1218 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1219 { 1220 PetscFunctionBegin; 1221 *flg = PETSC_TRUE; 1222 PetscFunctionReturn(0); 1223 } 1224 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1227 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1228 { 1229 PetscFunctionBegin; 1230 *flg = PETSC_FALSE; 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNCT__ 1235 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1236 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1237 { 1238 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1239 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1240 MatScalar *aa = a->a; 1241 1242 PetscFunctionBegin; 1243 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1244 PetscFunctionReturn(0); 1245 } 1246 1247 #undef __FUNCT__ 1248 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1249 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1250 { 1251 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1252 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1253 MatScalar *aa = a->a; 1254 1255 PetscFunctionBegin; 1256 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1257 PetscFunctionReturn(0); 1258 } 1259 1260 /* -------------------------------------------------------------------*/ 1261 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1262 MatGetRow_SeqSBAIJ, 1263 MatRestoreRow_SeqSBAIJ, 1264 MatMult_SeqSBAIJ_N, 1265 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1266 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1267 MatMultAdd_SeqSBAIJ_N, 1268 0, 1269 0, 1270 0, 1271 /*10*/ 0, 1272 0, 1273 MatCholeskyFactor_SeqSBAIJ, 1274 MatSOR_SeqSBAIJ, 1275 MatTranspose_SeqSBAIJ, 1276 /*15*/ MatGetInfo_SeqSBAIJ, 1277 MatEqual_SeqSBAIJ, 1278 MatGetDiagonal_SeqSBAIJ, 1279 MatDiagonalScale_SeqSBAIJ, 1280 MatNorm_SeqSBAIJ, 1281 /*20*/ 0, 1282 MatAssemblyEnd_SeqSBAIJ, 1283 MatSetOption_SeqSBAIJ, 1284 MatZeroEntries_SeqSBAIJ, 1285 /*24*/ 0, 1286 0, 1287 0, 1288 0, 1289 0, 1290 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1291 0, 1292 0, 1293 MatGetArray_SeqSBAIJ, 1294 MatRestoreArray_SeqSBAIJ, 1295 /*34*/ MatDuplicate_SeqSBAIJ, 1296 0, 1297 0, 1298 0, 1299 MatICCFactor_SeqSBAIJ, 1300 /*39*/ MatAXPY_SeqSBAIJ, 1301 MatGetSubMatrices_SeqSBAIJ, 1302 MatIncreaseOverlap_SeqSBAIJ, 1303 MatGetValues_SeqSBAIJ, 1304 MatCopy_SeqSBAIJ, 1305 /*44*/ 0, 1306 MatScale_SeqSBAIJ, 1307 0, 1308 0, 1309 0, 1310 /*49*/ 0, 1311 MatGetRowIJ_SeqSBAIJ, 1312 MatRestoreRowIJ_SeqSBAIJ, 1313 0, 1314 0, 1315 /*54*/ 0, 1316 0, 1317 0, 1318 0, 1319 MatSetValuesBlocked_SeqSBAIJ, 1320 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1321 0, 1322 0, 1323 0, 1324 0, 1325 /*64*/ 0, 1326 0, 1327 0, 1328 0, 1329 0, 1330 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1331 0, 1332 0, 1333 0, 1334 0, 1335 /*74*/ 0, 1336 0, 1337 0, 1338 0, 1339 0, 1340 /*79*/ 0, 1341 0, 1342 0, 1343 #if !defined(PETSC_USE_COMPLEX) 1344 MatGetInertia_SeqSBAIJ, 1345 #else 1346 0, 1347 #endif 1348 MatLoad_SeqSBAIJ, 1349 /*84*/ MatIsSymmetric_SeqSBAIJ, 1350 MatIsHermitian_SeqSBAIJ, 1351 MatIsStructurallySymmetric_SeqSBAIJ, 1352 0, 1353 0, 1354 /*89*/ 0, 1355 0, 1356 0, 1357 0, 1358 0, 1359 /*94*/ 0, 1360 0, 1361 0, 1362 0, 1363 0, 1364 /*99*/ 0, 1365 0, 1366 0, 1367 0, 1368 0, 1369 /*104*/0, 1370 MatRealPart_SeqSBAIJ, 1371 MatImaginaryPart_SeqSBAIJ, 1372 MatGetRowUpperTriangular_SeqSBAIJ, 1373 MatRestoreRowUpperTriangular_SeqSBAIJ, 1374 /*109*/0, 1375 0, 1376 0, 1377 0, 1378 MatMissingDiagonal_SeqSBAIJ 1379 }; 1380 1381 EXTERN_C_BEGIN 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1384 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1385 { 1386 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1387 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1388 PetscErrorCode ierr; 1389 1390 PetscFunctionBegin; 1391 if (aij->nonew != 1) { 1392 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1393 } 1394 1395 /* allocate space for values if not already there */ 1396 if (!aij->saved_values) { 1397 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1398 } 1399 1400 /* copy values over */ 1401 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1402 PetscFunctionReturn(0); 1403 } 1404 EXTERN_C_END 1405 1406 EXTERN_C_BEGIN 1407 #undef __FUNCT__ 1408 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1409 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1410 { 1411 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1412 PetscErrorCode ierr; 1413 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1414 1415 PetscFunctionBegin; 1416 if (aij->nonew != 1) { 1417 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1418 } 1419 if (!aij->saved_values) { 1420 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1421 } 1422 1423 /* copy values over */ 1424 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1425 PetscFunctionReturn(0); 1426 } 1427 EXTERN_C_END 1428 1429 EXTERN_C_BEGIN 1430 #undef __FUNCT__ 1431 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1432 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1433 { 1434 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1435 PetscErrorCode ierr; 1436 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1437 PetscTruth skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1438 1439 PetscFunctionBegin; 1440 B->preallocated = PETSC_TRUE; 1441 if (bs < 0) { 1442 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1443 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1444 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1445 bs = PetscAbs(bs); 1446 } 1447 if (nnz && newbs != bs) { 1448 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1449 } 1450 bs = newbs; 1451 1452 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1453 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1454 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1455 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1456 1457 mbs = B->rmap->N/bs; 1458 bs2 = bs*bs; 1459 1460 if (mbs*bs != B->rmap->N) { 1461 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1462 } 1463 1464 if (nz == MAT_SKIP_ALLOCATION) { 1465 skipallocation = PETSC_TRUE; 1466 nz = 0; 1467 } 1468 1469 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1470 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1471 if (nnz) { 1472 for (i=0; i<mbs; i++) { 1473 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1474 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); 1475 } 1476 } 1477 1478 B->ops->mult = MatMult_SeqSBAIJ_N; 1479 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1480 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1481 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1482 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1483 if (!flg) { 1484 switch (bs) { 1485 case 1: 1486 B->ops->mult = MatMult_SeqSBAIJ_1; 1487 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1488 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1489 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1490 break; 1491 case 2: 1492 B->ops->mult = MatMult_SeqSBAIJ_2; 1493 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1494 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1495 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1496 break; 1497 case 3: 1498 B->ops->mult = MatMult_SeqSBAIJ_3; 1499 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1500 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1501 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1502 break; 1503 case 4: 1504 B->ops->mult = MatMult_SeqSBAIJ_4; 1505 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1506 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1507 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1508 break; 1509 case 5: 1510 B->ops->mult = MatMult_SeqSBAIJ_5; 1511 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1512 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1513 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1514 break; 1515 case 6: 1516 B->ops->mult = MatMult_SeqSBAIJ_6; 1517 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1518 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1519 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1520 break; 1521 case 7: 1522 B->ops->mult = MatMult_SeqSBAIJ_7; 1523 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1524 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1525 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1526 break; 1527 } 1528 } 1529 1530 b->mbs = mbs; 1531 b->nbs = mbs; 1532 if (!skipallocation) { 1533 if (!b->imax) { 1534 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1535 b->free_imax_ilen = PETSC_TRUE; 1536 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1537 } 1538 if (!nnz) { 1539 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1540 else if (nz <= 0) nz = 1; 1541 for (i=0; i<mbs; i++) { 1542 b->imax[i] = nz; 1543 } 1544 nz = nz*mbs; /* total nz */ 1545 } else { 1546 nz = 0; 1547 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1548 } 1549 /* b->ilen will count nonzeros in each block row so far. */ 1550 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1551 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1552 1553 /* allocate the matrix space */ 1554 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1555 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1556 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1557 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1558 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1559 b->singlemalloc = PETSC_TRUE; 1560 1561 /* pointer to beginning of each row */ 1562 b->i[0] = 0; 1563 for (i=1; i<mbs+1; i++) { 1564 b->i[i] = b->i[i-1] + b->imax[i-1]; 1565 } 1566 b->free_a = PETSC_TRUE; 1567 b->free_ij = PETSC_TRUE; 1568 } else { 1569 b->free_a = PETSC_FALSE; 1570 b->free_ij = PETSC_FALSE; 1571 } 1572 1573 B->rmap->bs = bs; 1574 b->bs2 = bs2; 1575 b->nz = 0; 1576 b->maxnz = nz*bs2; 1577 1578 b->inew = 0; 1579 b->jnew = 0; 1580 b->anew = 0; 1581 b->a2anew = 0; 1582 b->permute = PETSC_FALSE; 1583 PetscFunctionReturn(0); 1584 } 1585 EXTERN_C_END 1586 1587 #undef __FUNCT__ 1588 #define __FUNCT__ "MatSeqBAIJSetNumericFactorization" 1589 /* 1590 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1591 */ 1592 PetscErrorCode MatSeqSBAIJSetNumericFactorization(Mat B,PetscTruth natural) 1593 { 1594 PetscErrorCode ierr; 1595 PetscTruth flg = PETSC_FALSE; 1596 PetscInt bs = B->rmap->bs; 1597 1598 PetscFunctionBegin; 1599 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1600 if (flg) bs = 8; 1601 1602 if (!natural) { 1603 switch (bs) { 1604 case 1: 1605 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1606 break; 1607 case 2: 1608 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1609 break; 1610 case 3: 1611 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1612 break; 1613 case 4: 1614 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1615 break; 1616 case 5: 1617 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1618 break; 1619 case 6: 1620 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1621 break; 1622 case 7: 1623 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1624 break; 1625 default: 1626 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1627 break; 1628 } 1629 } else { 1630 switch (bs) { 1631 case 1: 1632 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering; 1633 break; 1634 case 2: 1635 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1636 break; 1637 case 3: 1638 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1639 break; 1640 case 4: 1641 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1642 break; 1643 case 5: 1644 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1645 break; 1646 case 6: 1647 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1648 break; 1649 case 7: 1650 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1651 break; 1652 default: 1653 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1654 break; 1655 } 1656 } 1657 PetscFunctionReturn(0); 1658 } 1659 1660 EXTERN_C_BEGIN 1661 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1662 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1663 EXTERN_C_END 1664 1665 1666 EXTERN_C_BEGIN 1667 #undef __FUNCT__ 1668 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1669 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1670 { 1671 PetscInt n = A->rmap->n; 1672 PetscErrorCode ierr; 1673 1674 PetscFunctionBegin; 1675 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1676 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1677 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1678 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1679 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1680 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1681 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1682 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1683 (*B)->factor = ftype; 1684 PetscFunctionReturn(0); 1685 } 1686 EXTERN_C_END 1687 1688 EXTERN_C_BEGIN 1689 #undef __FUNCT__ 1690 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1691 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1692 { 1693 PetscFunctionBegin; 1694 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1695 *flg = PETSC_TRUE; 1696 } else { 1697 *flg = PETSC_FALSE; 1698 } 1699 PetscFunctionReturn(0); 1700 } 1701 EXTERN_C_END 1702 1703 EXTERN_C_BEGIN 1704 #if defined(PETSC_HAVE_MUMPS) 1705 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1706 #endif 1707 #if defined(PETSC_HAVE_SPOOLES) 1708 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1709 #endif 1710 #if defined(PETSC_HAVE_PASTIX) 1711 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1712 #endif 1713 EXTERN_C_END 1714 1715 /*MC 1716 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1717 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1718 1719 Options Database Keys: 1720 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1721 1722 Level: beginner 1723 1724 .seealso: MatCreateSeqSBAIJ 1725 M*/ 1726 1727 EXTERN_C_BEGIN 1728 #undef __FUNCT__ 1729 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1730 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1731 { 1732 Mat_SeqSBAIJ *b; 1733 PetscErrorCode ierr; 1734 PetscMPIInt size; 1735 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1736 1737 PetscFunctionBegin; 1738 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1739 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1740 1741 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1742 B->data = (void*)b; 1743 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1744 B->ops->destroy = MatDestroy_SeqSBAIJ; 1745 B->ops->view = MatView_SeqSBAIJ; 1746 B->mapping = 0; 1747 b->row = 0; 1748 b->icol = 0; 1749 b->reallocs = 0; 1750 b->saved_values = 0; 1751 b->inode.limit = 5; 1752 b->inode.max_limit = 5; 1753 1754 b->roworiented = PETSC_TRUE; 1755 b->nonew = 0; 1756 b->diag = 0; 1757 b->solve_work = 0; 1758 b->mult_work = 0; 1759 B->spptr = 0; 1760 b->keepnonzeropattern = PETSC_FALSE; 1761 b->xtoy = 0; 1762 b->XtoY = 0; 1763 1764 b->inew = 0; 1765 b->jnew = 0; 1766 b->anew = 0; 1767 b->a2anew = 0; 1768 b->permute = PETSC_FALSE; 1769 1770 b->ignore_ltriangular = PETSC_FALSE; 1771 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1772 1773 b->getrow_utriangular = PETSC_FALSE; 1774 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1775 1776 #if defined(PETSC_HAVE_PASTIX) 1777 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1778 "MatGetFactor_seqsbaij_pastix", 1779 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1780 #endif 1781 #if defined(PETSC_HAVE_SPOOLES) 1782 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1783 "MatGetFactor_seqsbaij_spooles", 1784 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1785 #endif 1786 #if defined(PETSC_HAVE_MUMPS) 1787 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1788 "MatGetFactor_seqsbaij_mumps", 1789 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1790 #endif 1791 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1792 "MatGetFactorAvailable_seqsbaij_petsc", 1793 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1794 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1795 "MatGetFactor_seqsbaij_petsc", 1796 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1797 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1798 "MatStoreValues_SeqSBAIJ", 1799 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1800 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1801 "MatRetrieveValues_SeqSBAIJ", 1802 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1803 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1804 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1805 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1806 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1807 "MatConvert_SeqSBAIJ_SeqAIJ", 1808 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1809 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1810 "MatConvert_SeqSBAIJ_SeqBAIJ", 1811 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1812 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1813 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1814 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1815 1816 B->symmetric = PETSC_TRUE; 1817 B->structurally_symmetric = PETSC_TRUE; 1818 B->symmetric_set = PETSC_TRUE; 1819 B->structurally_symmetric_set = PETSC_TRUE; 1820 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1821 1822 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1823 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1824 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1825 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1826 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1827 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); 1828 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1829 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1830 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1831 1832 PetscFunctionReturn(0); 1833 } 1834 EXTERN_C_END 1835 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1838 /*@C 1839 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1840 compressed row) format. For good matrix assembly performance the 1841 user should preallocate the matrix storage by setting the parameter nz 1842 (or the array nnz). By setting these parameters accurately, performance 1843 during matrix assembly can be increased by more than a factor of 50. 1844 1845 Collective on Mat 1846 1847 Input Parameters: 1848 + A - the symmetric matrix 1849 . bs - size of block 1850 . nz - number of block nonzeros per block row (same for all rows) 1851 - nnz - array containing the number of block nonzeros in the upper triangular plus 1852 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1853 1854 Options Database Keys: 1855 . -mat_no_unroll - uses code that does not unroll the loops in the 1856 block calculations (much slower) 1857 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1858 1859 Level: intermediate 1860 1861 Notes: 1862 Specify the preallocated storage with either nz or nnz (not both). 1863 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1864 allocation. For additional details, see the users manual chapter on 1865 matrices. 1866 1867 You can call MatGetInfo() to get information on how effective the preallocation was; 1868 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1869 You can also run with the option -info and look for messages with the string 1870 malloc in them to see if additional memory allocation was needed. 1871 1872 If the nnz parameter is given then the nz parameter is ignored 1873 1874 1875 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1876 @*/ 1877 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1878 { 1879 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1880 1881 PetscFunctionBegin; 1882 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1883 if (f) { 1884 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1885 } 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "MatCreateSeqSBAIJ" 1891 /*@C 1892 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1893 compressed row) format. For good matrix assembly performance the 1894 user should preallocate the matrix storage by setting the parameter nz 1895 (or the array nnz). By setting these parameters accurately, performance 1896 during matrix assembly can be increased by more than a factor of 50. 1897 1898 Collective on MPI_Comm 1899 1900 Input Parameters: 1901 + comm - MPI communicator, set to PETSC_COMM_SELF 1902 . bs - size of block 1903 . m - number of rows, or number of columns 1904 . nz - number of block nonzeros per block row (same for all rows) 1905 - nnz - array containing the number of block nonzeros in the upper triangular plus 1906 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1907 1908 Output Parameter: 1909 . A - the symmetric matrix 1910 1911 Options Database Keys: 1912 . -mat_no_unroll - uses code that does not unroll the loops in the 1913 block calculations (much slower) 1914 . -mat_block_size - size of the blocks to use 1915 1916 Level: intermediate 1917 1918 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1919 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1920 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1921 1922 Notes: 1923 The number of rows and columns must be divisible by blocksize. 1924 This matrix type does not support complex Hermitian operation. 1925 1926 Specify the preallocated storage with either nz or nnz (not both). 1927 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1928 allocation. For additional details, see the users manual chapter on 1929 matrices. 1930 1931 If the nnz parameter is given then the nz parameter is ignored 1932 1933 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1934 @*/ 1935 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1936 { 1937 PetscErrorCode ierr; 1938 1939 PetscFunctionBegin; 1940 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1941 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1942 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1943 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 #undef __FUNCT__ 1948 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1949 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1950 { 1951 Mat C; 1952 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1953 PetscErrorCode ierr; 1954 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1955 1956 PetscFunctionBegin; 1957 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1958 1959 *B = 0; 1960 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1961 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1962 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1963 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1964 c = (Mat_SeqSBAIJ*)C->data; 1965 1966 C->preallocated = PETSC_TRUE; 1967 C->factor = A->factor; 1968 c->row = 0; 1969 c->icol = 0; 1970 c->saved_values = 0; 1971 c->keepnonzeropattern = a->keepnonzeropattern; 1972 C->assembled = PETSC_TRUE; 1973 1974 ierr = PetscLayoutCopy(A->rmap,&C->rmap);CHKERRQ(ierr); 1975 ierr = PetscLayoutCopy(A->cmap,&C->cmap);CHKERRQ(ierr); 1976 c->bs2 = a->bs2; 1977 c->mbs = a->mbs; 1978 c->nbs = a->nbs; 1979 1980 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 1981 c->imax = a->imax; 1982 c->ilen = a->ilen; 1983 c->free_imax_ilen = PETSC_FALSE; 1984 } else { 1985 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1986 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1987 for (i=0; i<mbs; i++) { 1988 c->imax[i] = a->imax[i]; 1989 c->ilen[i] = a->ilen[i]; 1990 } 1991 c->free_imax_ilen = PETSC_TRUE; 1992 } 1993 1994 /* allocate the matrix space */ 1995 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 1996 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 1997 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1998 c->singlemalloc = PETSC_FALSE; 1999 c->free_ij = PETSC_FALSE; 2000 c->parent = A; 2001 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2002 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2003 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2004 } else { 2005 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2006 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2007 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2008 c->singlemalloc = PETSC_TRUE; 2009 c->free_ij = PETSC_TRUE; 2010 } 2011 if (mbs > 0) { 2012 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2013 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2014 } 2015 if (cpvalues == MAT_COPY_VALUES) { 2016 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2017 } else { 2018 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2019 } 2020 if (a->jshort) { 2021 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2022 c->jshort = a->jshort; 2023 c->free_jshort = PETSC_FALSE; 2024 } else { 2025 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2026 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2027 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2028 c->free_jshort = PETSC_TRUE; 2029 } 2030 } 2031 } 2032 2033 c->roworiented = a->roworiented; 2034 c->nonew = a->nonew; 2035 2036 if (a->diag) { 2037 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2038 c->diag = a->diag; 2039 c->free_diag = PETSC_FALSE; 2040 } else { 2041 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2042 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2043 for (i=0; i<mbs; i++) { 2044 c->diag[i] = a->diag[i]; 2045 } 2046 c->free_diag = PETSC_TRUE; 2047 } 2048 } else c->diag = 0; 2049 c->nz = a->nz; 2050 c->maxnz = a->maxnz; 2051 c->solve_work = 0; 2052 c->mult_work = 0; 2053 c->free_a = PETSC_TRUE; 2054 *B = C; 2055 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2056 PetscFunctionReturn(0); 2057 } 2058 2059 #undef __FUNCT__ 2060 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2061 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 2062 { 2063 Mat_SeqSBAIJ *a; 2064 Mat B; 2065 PetscErrorCode ierr; 2066 int fd; 2067 PetscMPIInt size; 2068 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2069 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2070 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2071 PetscInt *masked,nmask,tmp,bs2,ishift; 2072 PetscScalar *aa; 2073 MPI_Comm comm = ((PetscObject)viewer)->comm; 2074 2075 PetscFunctionBegin; 2076 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2077 bs2 = bs*bs; 2078 2079 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2080 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 2081 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2082 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2083 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2084 M = header[1]; N = header[2]; nz = header[3]; 2085 2086 if (header[3] < 0) { 2087 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2088 } 2089 2090 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 2091 2092 /* 2093 This code adds extra rows to make sure the number of rows is 2094 divisible by the blocksize 2095 */ 2096 mbs = M/bs; 2097 extra_rows = bs - M + bs*(mbs); 2098 if (extra_rows == bs) extra_rows = 0; 2099 else mbs++; 2100 if (extra_rows) { 2101 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2102 } 2103 2104 /* read in row lengths */ 2105 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2106 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2107 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2108 2109 /* read in column indices */ 2110 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2111 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2112 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2113 2114 /* loop over row lengths determining block row lengths */ 2115 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2116 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2117 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 2118 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2119 masked = mask + mbs; 2120 rowcount = 0; nzcount = 0; 2121 for (i=0; i<mbs; i++) { 2122 nmask = 0; 2123 for (j=0; j<bs; j++) { 2124 kmax = rowlengths[rowcount]; 2125 for (k=0; k<kmax; k++) { 2126 tmp = jj[nzcount++]/bs; /* block col. index */ 2127 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2128 } 2129 rowcount++; 2130 } 2131 s_browlengths[i] += nmask; 2132 2133 /* zero out the mask elements we set */ 2134 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2135 } 2136 2137 /* create our matrix */ 2138 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2139 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2140 ierr = MatSetType(B,type);CHKERRQ(ierr); 2141 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2142 a = (Mat_SeqSBAIJ*)B->data; 2143 2144 /* set matrix "i" values */ 2145 a->i[0] = 0; 2146 for (i=1; i<= mbs; i++) { 2147 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2148 a->ilen[i-1] = s_browlengths[i-1]; 2149 } 2150 a->nz = a->i[mbs]; 2151 2152 /* read in nonzero values */ 2153 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2154 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2155 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2156 2157 /* set "a" and "j" values into matrix */ 2158 nzcount = 0; jcount = 0; 2159 for (i=0; i<mbs; i++) { 2160 nzcountb = nzcount; 2161 nmask = 0; 2162 for (j=0; j<bs; j++) { 2163 kmax = rowlengths[i*bs+j]; 2164 for (k=0; k<kmax; k++) { 2165 tmp = jj[nzcount++]/bs; /* block col. index */ 2166 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2167 } 2168 } 2169 /* sort the masked values */ 2170 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2171 2172 /* set "j" values into matrix */ 2173 maskcount = 1; 2174 for (j=0; j<nmask; j++) { 2175 a->j[jcount++] = masked[j]; 2176 mask[masked[j]] = maskcount++; 2177 } 2178 2179 /* set "a" values into matrix */ 2180 ishift = bs2*a->i[i]; 2181 for (j=0; j<bs; j++) { 2182 kmax = rowlengths[i*bs+j]; 2183 for (k=0; k<kmax; k++) { 2184 tmp = jj[nzcountb]/bs ; /* block col. index */ 2185 if (tmp >= i){ 2186 block = mask[tmp] - 1; 2187 point = jj[nzcountb] - bs*tmp; 2188 idx = ishift + bs2*block + j + bs*point; 2189 a->a[idx] = aa[nzcountb]; 2190 } 2191 nzcountb++; 2192 } 2193 } 2194 /* zero out the mask elements we set */ 2195 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2196 } 2197 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2198 2199 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2200 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2201 ierr = PetscFree(aa);CHKERRQ(ierr); 2202 ierr = PetscFree(jj);CHKERRQ(ierr); 2203 ierr = PetscFree(mask);CHKERRQ(ierr); 2204 2205 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2206 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2207 ierr = MatView_Private(B);CHKERRQ(ierr); 2208 *A = B; 2209 PetscFunctionReturn(0); 2210 } 2211 2212 #undef __FUNCT__ 2213 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2214 /*@ 2215 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2216 (upper triangular entries in CSR format) provided by the user. 2217 2218 Collective on MPI_Comm 2219 2220 Input Parameters: 2221 + comm - must be an MPI communicator of size 1 2222 . bs - size of block 2223 . m - number of rows 2224 . n - number of columns 2225 . i - row indices 2226 . j - column indices 2227 - a - matrix values 2228 2229 Output Parameter: 2230 . mat - the matrix 2231 2232 Level: intermediate 2233 2234 Notes: 2235 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2236 once the matrix is destroyed 2237 2238 You cannot set new nonzero locations into this matrix, that will generate an error. 2239 2240 The i and j indices are 0 based 2241 2242 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2243 2244 @*/ 2245 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2246 { 2247 PetscErrorCode ierr; 2248 PetscInt ii; 2249 Mat_SeqSBAIJ *sbaij; 2250 2251 PetscFunctionBegin; 2252 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2253 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2254 2255 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2256 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2257 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2258 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2259 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2260 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2261 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2262 2263 sbaij->i = i; 2264 sbaij->j = j; 2265 sbaij->a = a; 2266 sbaij->singlemalloc = PETSC_FALSE; 2267 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2268 sbaij->free_a = PETSC_FALSE; 2269 sbaij->free_ij = PETSC_FALSE; 2270 2271 for (ii=0; ii<m; ii++) { 2272 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2273 #if defined(PETSC_USE_DEBUG) 2274 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]); 2275 #endif 2276 } 2277 #if defined(PETSC_USE_DEBUG) 2278 for (ii=0; ii<sbaij->i[m]; ii++) { 2279 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2280 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2281 } 2282 #endif 2283 2284 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2285 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2286 PetscFunctionReturn(0); 2287 } 2288 2289 2290 2291 2292 2293