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