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 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 /*114*/0, 1380 0, 1381 0, 1382 0, 1383 0, 1384 /*119*/0, 1385 0, 1386 0 1387 }; 1388 1389 EXTERN_C_BEGIN 1390 #undef __FUNCT__ 1391 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1392 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1393 { 1394 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1395 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1396 PetscErrorCode ierr; 1397 1398 PetscFunctionBegin; 1399 if (aij->nonew != 1) { 1400 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1401 } 1402 1403 /* allocate space for values if not already there */ 1404 if (!aij->saved_values) { 1405 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1406 } 1407 1408 /* copy values over */ 1409 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 EXTERN_C_END 1413 1414 EXTERN_C_BEGIN 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1417 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1418 { 1419 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1420 PetscErrorCode ierr; 1421 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1422 1423 PetscFunctionBegin; 1424 if (aij->nonew != 1) { 1425 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1426 } 1427 if (!aij->saved_values) { 1428 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1429 } 1430 1431 /* copy values over */ 1432 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1433 PetscFunctionReturn(0); 1434 } 1435 EXTERN_C_END 1436 1437 EXTERN_C_BEGIN 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1440 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1441 { 1442 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1443 PetscErrorCode ierr; 1444 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1445 PetscTruth skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1446 1447 PetscFunctionBegin; 1448 B->preallocated = PETSC_TRUE; 1449 if (bs < 0) { 1450 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1451 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1452 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1453 bs = PetscAbs(bs); 1454 } 1455 if (nnz && newbs != bs) { 1456 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1457 } 1458 bs = newbs; 1459 1460 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1461 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1462 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1463 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1464 1465 mbs = B->rmap->N/bs; 1466 bs2 = bs*bs; 1467 1468 if (mbs*bs != B->rmap->N) { 1469 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1470 } 1471 1472 if (nz == MAT_SKIP_ALLOCATION) { 1473 skipallocation = PETSC_TRUE; 1474 nz = 0; 1475 } 1476 1477 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1478 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1479 if (nnz) { 1480 for (i=0; i<mbs; i++) { 1481 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1482 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); 1483 } 1484 } 1485 1486 B->ops->mult = MatMult_SeqSBAIJ_N; 1487 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1488 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1489 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1490 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1491 if (!flg) { 1492 switch (bs) { 1493 case 1: 1494 B->ops->mult = MatMult_SeqSBAIJ_1; 1495 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1496 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1497 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1498 break; 1499 case 2: 1500 B->ops->mult = MatMult_SeqSBAIJ_2; 1501 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1502 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1503 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1504 break; 1505 case 3: 1506 B->ops->mult = MatMult_SeqSBAIJ_3; 1507 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1508 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1509 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1510 break; 1511 case 4: 1512 B->ops->mult = MatMult_SeqSBAIJ_4; 1513 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1514 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1515 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1516 break; 1517 case 5: 1518 B->ops->mult = MatMult_SeqSBAIJ_5; 1519 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1520 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1521 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1522 break; 1523 case 6: 1524 B->ops->mult = MatMult_SeqSBAIJ_6; 1525 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1526 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1527 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1528 break; 1529 case 7: 1530 B->ops->mult = MatMult_SeqSBAIJ_7; 1531 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1532 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1533 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1534 break; 1535 } 1536 } 1537 1538 b->mbs = mbs; 1539 b->nbs = mbs; 1540 if (!skipallocation) { 1541 if (!b->imax) { 1542 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1543 b->free_imax_ilen = PETSC_TRUE; 1544 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1545 } 1546 if (!nnz) { 1547 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1548 else if (nz <= 0) nz = 1; 1549 for (i=0; i<mbs; i++) { 1550 b->imax[i] = nz; 1551 } 1552 nz = nz*mbs; /* total nz */ 1553 } else { 1554 nz = 0; 1555 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1556 } 1557 /* b->ilen will count nonzeros in each block row so far. */ 1558 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1559 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1560 1561 /* allocate the matrix space */ 1562 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1563 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1564 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1565 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1566 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1567 b->singlemalloc = PETSC_TRUE; 1568 1569 /* pointer to beginning of each row */ 1570 b->i[0] = 0; 1571 for (i=1; i<mbs+1; i++) { 1572 b->i[i] = b->i[i-1] + b->imax[i-1]; 1573 } 1574 b->free_a = PETSC_TRUE; 1575 b->free_ij = PETSC_TRUE; 1576 } else { 1577 b->free_a = PETSC_FALSE; 1578 b->free_ij = PETSC_FALSE; 1579 } 1580 1581 B->rmap->bs = bs; 1582 b->bs2 = bs2; 1583 b->nz = 0; 1584 b->maxnz = nz*bs2; 1585 1586 b->inew = 0; 1587 b->jnew = 0; 1588 b->anew = 0; 1589 b->a2anew = 0; 1590 b->permute = PETSC_FALSE; 1591 PetscFunctionReturn(0); 1592 } 1593 EXTERN_C_END 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "MatSeqBAIJSetNumericFactorization" 1597 /* 1598 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1599 */ 1600 PetscErrorCode MatSeqSBAIJSetNumericFactorization(Mat B,PetscTruth natural) 1601 { 1602 PetscErrorCode ierr; 1603 PetscTruth flg = PETSC_FALSE; 1604 PetscInt bs = B->rmap->bs; 1605 1606 PetscFunctionBegin; 1607 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1608 if (flg) bs = 8; 1609 1610 if (!natural) { 1611 switch (bs) { 1612 case 1: 1613 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1614 break; 1615 case 2: 1616 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1617 break; 1618 case 3: 1619 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1620 break; 1621 case 4: 1622 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1623 break; 1624 case 5: 1625 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1626 break; 1627 case 6: 1628 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1629 break; 1630 case 7: 1631 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1632 break; 1633 default: 1634 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1635 break; 1636 } 1637 } else { 1638 switch (bs) { 1639 case 1: 1640 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering; 1641 break; 1642 case 2: 1643 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1644 break; 1645 case 3: 1646 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1647 break; 1648 case 4: 1649 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1650 break; 1651 case 5: 1652 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1653 break; 1654 case 6: 1655 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1656 break; 1657 case 7: 1658 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1659 break; 1660 default: 1661 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1662 break; 1663 } 1664 } 1665 PetscFunctionReturn(0); 1666 } 1667 1668 EXTERN_C_BEGIN 1669 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1670 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1671 EXTERN_C_END 1672 1673 1674 EXTERN_C_BEGIN 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1677 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1678 { 1679 PetscInt n = A->rmap->n; 1680 PetscErrorCode ierr; 1681 1682 PetscFunctionBegin; 1683 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1684 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1685 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1686 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1687 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1688 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1689 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1690 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1691 (*B)->factor = ftype; 1692 PetscFunctionReturn(0); 1693 } 1694 EXTERN_C_END 1695 1696 EXTERN_C_BEGIN 1697 #undef __FUNCT__ 1698 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1699 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1700 { 1701 PetscFunctionBegin; 1702 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1703 *flg = PETSC_TRUE; 1704 } else { 1705 *flg = PETSC_FALSE; 1706 } 1707 PetscFunctionReturn(0); 1708 } 1709 EXTERN_C_END 1710 1711 EXTERN_C_BEGIN 1712 #if defined(PETSC_HAVE_MUMPS) 1713 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1714 #endif 1715 #if defined(PETSC_HAVE_SPOOLES) 1716 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1717 #endif 1718 #if defined(PETSC_HAVE_PASTIX) 1719 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1720 #endif 1721 EXTERN_C_END 1722 1723 /*MC 1724 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1725 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1726 1727 Options Database Keys: 1728 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1729 1730 Level: beginner 1731 1732 .seealso: MatCreateSeqSBAIJ 1733 M*/ 1734 1735 EXTERN_C_BEGIN 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1738 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1739 { 1740 Mat_SeqSBAIJ *b; 1741 PetscErrorCode ierr; 1742 PetscMPIInt size; 1743 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1744 1745 PetscFunctionBegin; 1746 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1747 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1748 1749 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1750 B->data = (void*)b; 1751 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1752 B->ops->destroy = MatDestroy_SeqSBAIJ; 1753 B->ops->view = MatView_SeqSBAIJ; 1754 B->mapping = 0; 1755 b->row = 0; 1756 b->icol = 0; 1757 b->reallocs = 0; 1758 b->saved_values = 0; 1759 b->inode.limit = 5; 1760 b->inode.max_limit = 5; 1761 1762 b->roworiented = PETSC_TRUE; 1763 b->nonew = 0; 1764 b->diag = 0; 1765 b->solve_work = 0; 1766 b->mult_work = 0; 1767 B->spptr = 0; 1768 b->keepnonzeropattern = PETSC_FALSE; 1769 b->xtoy = 0; 1770 b->XtoY = 0; 1771 1772 b->inew = 0; 1773 b->jnew = 0; 1774 b->anew = 0; 1775 b->a2anew = 0; 1776 b->permute = PETSC_FALSE; 1777 1778 b->ignore_ltriangular = PETSC_FALSE; 1779 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1780 1781 b->getrow_utriangular = PETSC_FALSE; 1782 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1783 1784 #if defined(PETSC_HAVE_PASTIX) 1785 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1786 "MatGetFactor_seqsbaij_pastix", 1787 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1788 #endif 1789 #if defined(PETSC_HAVE_SPOOLES) 1790 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1791 "MatGetFactor_seqsbaij_spooles", 1792 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1793 #endif 1794 #if defined(PETSC_HAVE_MUMPS) 1795 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1796 "MatGetFactor_seqsbaij_mumps", 1797 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1798 #endif 1799 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1800 "MatGetFactorAvailable_seqsbaij_petsc", 1801 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1802 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1803 "MatGetFactor_seqsbaij_petsc", 1804 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1805 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1806 "MatStoreValues_SeqSBAIJ", 1807 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1808 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1809 "MatRetrieveValues_SeqSBAIJ", 1810 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1811 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1812 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1813 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1814 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1815 "MatConvert_SeqSBAIJ_SeqAIJ", 1816 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1817 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1818 "MatConvert_SeqSBAIJ_SeqBAIJ", 1819 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1821 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1822 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1823 1824 B->symmetric = PETSC_TRUE; 1825 B->structurally_symmetric = PETSC_TRUE; 1826 B->symmetric_set = PETSC_TRUE; 1827 B->structurally_symmetric_set = PETSC_TRUE; 1828 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1829 1830 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1831 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1832 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1833 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1834 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1835 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); 1836 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1837 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1838 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1839 1840 PetscFunctionReturn(0); 1841 } 1842 EXTERN_C_END 1843 1844 #undef __FUNCT__ 1845 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1846 /*@C 1847 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1848 compressed row) format. For good matrix assembly performance the 1849 user should preallocate the matrix storage by setting the parameter nz 1850 (or the array nnz). By setting these parameters accurately, performance 1851 during matrix assembly can be increased by more than a factor of 50. 1852 1853 Collective on Mat 1854 1855 Input Parameters: 1856 + A - the symmetric matrix 1857 . bs - size of block 1858 . nz - number of block nonzeros per block row (same for all rows) 1859 - nnz - array containing the number of block nonzeros in the upper triangular plus 1860 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1861 1862 Options Database Keys: 1863 . -mat_no_unroll - uses code that does not unroll the loops in the 1864 block calculations (much slower) 1865 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1866 1867 Level: intermediate 1868 1869 Notes: 1870 Specify the preallocated storage with either nz or nnz (not both). 1871 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1872 allocation. For additional details, see the users manual chapter on 1873 matrices. 1874 1875 You can call MatGetInfo() to get information on how effective the preallocation was; 1876 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1877 You can also run with the option -info and look for messages with the string 1878 malloc in them to see if additional memory allocation was needed. 1879 1880 If the nnz parameter is given then the nz parameter is ignored 1881 1882 1883 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1884 @*/ 1885 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1886 { 1887 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1888 1889 PetscFunctionBegin; 1890 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1891 if (f) { 1892 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1893 } 1894 PetscFunctionReturn(0); 1895 } 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatCreateSeqSBAIJ" 1899 /*@C 1900 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1901 compressed row) format. For good matrix assembly performance the 1902 user should preallocate the matrix storage by setting the parameter nz 1903 (or the array nnz). By setting these parameters accurately, performance 1904 during matrix assembly can be increased by more than a factor of 50. 1905 1906 Collective on MPI_Comm 1907 1908 Input Parameters: 1909 + comm - MPI communicator, set to PETSC_COMM_SELF 1910 . bs - size of block 1911 . m - number of rows, or number of columns 1912 . nz - number of block nonzeros per block row (same for all rows) 1913 - nnz - array containing the number of block nonzeros in the upper triangular plus 1914 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1915 1916 Output Parameter: 1917 . A - the symmetric matrix 1918 1919 Options Database Keys: 1920 . -mat_no_unroll - uses code that does not unroll the loops in the 1921 block calculations (much slower) 1922 . -mat_block_size - size of the blocks to use 1923 1924 Level: intermediate 1925 1926 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1927 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1928 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1929 1930 Notes: 1931 The number of rows and columns must be divisible by blocksize. 1932 This matrix type does not support complex Hermitian operation. 1933 1934 Specify the preallocated storage with either nz or nnz (not both). 1935 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1936 allocation. For additional details, see the users manual chapter on 1937 matrices. 1938 1939 If the nnz parameter is given then the nz parameter is ignored 1940 1941 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1942 @*/ 1943 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1944 { 1945 PetscErrorCode ierr; 1946 1947 PetscFunctionBegin; 1948 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1949 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1950 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1951 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1952 PetscFunctionReturn(0); 1953 } 1954 1955 #undef __FUNCT__ 1956 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1957 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1958 { 1959 Mat C; 1960 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1961 PetscErrorCode ierr; 1962 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1963 1964 PetscFunctionBegin; 1965 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1966 1967 *B = 0; 1968 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1969 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1970 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1971 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1972 c = (Mat_SeqSBAIJ*)C->data; 1973 1974 C->preallocated = PETSC_TRUE; 1975 C->factor = A->factor; 1976 c->row = 0; 1977 c->icol = 0; 1978 c->saved_values = 0; 1979 c->keepnonzeropattern = a->keepnonzeropattern; 1980 C->assembled = PETSC_TRUE; 1981 1982 ierr = PetscLayoutCopy(A->rmap,&C->rmap);CHKERRQ(ierr); 1983 ierr = PetscLayoutCopy(A->cmap,&C->cmap);CHKERRQ(ierr); 1984 c->bs2 = a->bs2; 1985 c->mbs = a->mbs; 1986 c->nbs = a->nbs; 1987 1988 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 1989 c->imax = a->imax; 1990 c->ilen = a->ilen; 1991 c->free_imax_ilen = PETSC_FALSE; 1992 } else { 1993 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1994 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1995 for (i=0; i<mbs; i++) { 1996 c->imax[i] = a->imax[i]; 1997 c->ilen[i] = a->ilen[i]; 1998 } 1999 c->free_imax_ilen = PETSC_TRUE; 2000 } 2001 2002 /* allocate the matrix space */ 2003 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2004 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2005 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2006 c->singlemalloc = PETSC_FALSE; 2007 c->free_ij = PETSC_FALSE; 2008 c->parent = A; 2009 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2010 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2011 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2012 } else { 2013 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2014 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2015 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2016 c->singlemalloc = PETSC_TRUE; 2017 c->free_ij = PETSC_TRUE; 2018 } 2019 if (mbs > 0) { 2020 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2021 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2022 } 2023 if (cpvalues == MAT_COPY_VALUES) { 2024 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2025 } else { 2026 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2027 } 2028 if (a->jshort) { 2029 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2030 c->jshort = a->jshort; 2031 c->free_jshort = PETSC_FALSE; 2032 } else { 2033 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2034 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2035 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2036 c->free_jshort = PETSC_TRUE; 2037 } 2038 } 2039 } 2040 2041 c->roworiented = a->roworiented; 2042 c->nonew = a->nonew; 2043 2044 if (a->diag) { 2045 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2046 c->diag = a->diag; 2047 c->free_diag = PETSC_FALSE; 2048 } else { 2049 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2050 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2051 for (i=0; i<mbs; i++) { 2052 c->diag[i] = a->diag[i]; 2053 } 2054 c->free_diag = PETSC_TRUE; 2055 } 2056 } else c->diag = 0; 2057 c->nz = a->nz; 2058 c->maxnz = a->maxnz; 2059 c->solve_work = 0; 2060 c->mult_work = 0; 2061 c->free_a = PETSC_TRUE; 2062 *B = C; 2063 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 #undef __FUNCT__ 2068 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2069 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 2070 { 2071 Mat_SeqSBAIJ *a; 2072 Mat B; 2073 PetscErrorCode ierr; 2074 int fd; 2075 PetscMPIInt size; 2076 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2077 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2078 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2079 PetscInt *masked,nmask,tmp,bs2,ishift; 2080 PetscScalar *aa; 2081 MPI_Comm comm = ((PetscObject)viewer)->comm; 2082 2083 PetscFunctionBegin; 2084 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2085 bs2 = bs*bs; 2086 2087 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2088 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 2089 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2090 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2091 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2092 M = header[1]; N = header[2]; nz = header[3]; 2093 2094 if (header[3] < 0) { 2095 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2096 } 2097 2098 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 2099 2100 /* 2101 This code adds extra rows to make sure the number of rows is 2102 divisible by the blocksize 2103 */ 2104 mbs = M/bs; 2105 extra_rows = bs - M + bs*(mbs); 2106 if (extra_rows == bs) extra_rows = 0; 2107 else mbs++; 2108 if (extra_rows) { 2109 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2110 } 2111 2112 /* read in row lengths */ 2113 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2114 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2115 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2116 2117 /* read in column indices */ 2118 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2119 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2120 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2121 2122 /* loop over row lengths determining block row lengths */ 2123 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2124 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2125 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 2126 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2127 masked = mask + mbs; 2128 rowcount = 0; nzcount = 0; 2129 for (i=0; i<mbs; i++) { 2130 nmask = 0; 2131 for (j=0; j<bs; j++) { 2132 kmax = rowlengths[rowcount]; 2133 for (k=0; k<kmax; k++) { 2134 tmp = jj[nzcount++]/bs; /* block col. index */ 2135 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2136 } 2137 rowcount++; 2138 } 2139 s_browlengths[i] += nmask; 2140 2141 /* zero out the mask elements we set */ 2142 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2143 } 2144 2145 /* create our matrix */ 2146 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2147 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2148 ierr = MatSetType(B,type);CHKERRQ(ierr); 2149 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2150 a = (Mat_SeqSBAIJ*)B->data; 2151 2152 /* set matrix "i" values */ 2153 a->i[0] = 0; 2154 for (i=1; i<= mbs; i++) { 2155 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2156 a->ilen[i-1] = s_browlengths[i-1]; 2157 } 2158 a->nz = a->i[mbs]; 2159 2160 /* read in nonzero values */ 2161 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2162 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2163 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2164 2165 /* set "a" and "j" values into matrix */ 2166 nzcount = 0; jcount = 0; 2167 for (i=0; i<mbs; i++) { 2168 nzcountb = nzcount; 2169 nmask = 0; 2170 for (j=0; j<bs; j++) { 2171 kmax = rowlengths[i*bs+j]; 2172 for (k=0; k<kmax; k++) { 2173 tmp = jj[nzcount++]/bs; /* block col. index */ 2174 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2175 } 2176 } 2177 /* sort the masked values */ 2178 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2179 2180 /* set "j" values into matrix */ 2181 maskcount = 1; 2182 for (j=0; j<nmask; j++) { 2183 a->j[jcount++] = masked[j]; 2184 mask[masked[j]] = maskcount++; 2185 } 2186 2187 /* set "a" values into matrix */ 2188 ishift = bs2*a->i[i]; 2189 for (j=0; j<bs; j++) { 2190 kmax = rowlengths[i*bs+j]; 2191 for (k=0; k<kmax; k++) { 2192 tmp = jj[nzcountb]/bs ; /* block col. index */ 2193 if (tmp >= i){ 2194 block = mask[tmp] - 1; 2195 point = jj[nzcountb] - bs*tmp; 2196 idx = ishift + bs2*block + j + bs*point; 2197 a->a[idx] = aa[nzcountb]; 2198 } 2199 nzcountb++; 2200 } 2201 } 2202 /* zero out the mask elements we set */ 2203 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2204 } 2205 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2206 2207 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2208 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2209 ierr = PetscFree(aa);CHKERRQ(ierr); 2210 ierr = PetscFree(jj);CHKERRQ(ierr); 2211 ierr = PetscFree(mask);CHKERRQ(ierr); 2212 2213 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2214 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2215 ierr = MatView_Private(B);CHKERRQ(ierr); 2216 *A = B; 2217 PetscFunctionReturn(0); 2218 } 2219 2220 #undef __FUNCT__ 2221 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2222 /*@ 2223 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2224 (upper triangular entries in CSR format) provided by the user. 2225 2226 Collective on MPI_Comm 2227 2228 Input Parameters: 2229 + comm - must be an MPI communicator of size 1 2230 . bs - size of block 2231 . m - number of rows 2232 . n - number of columns 2233 . i - row indices 2234 . j - column indices 2235 - a - matrix values 2236 2237 Output Parameter: 2238 . mat - the matrix 2239 2240 Level: intermediate 2241 2242 Notes: 2243 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2244 once the matrix is destroyed 2245 2246 You cannot set new nonzero locations into this matrix, that will generate an error. 2247 2248 The i and j indices are 0 based 2249 2250 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2251 2252 @*/ 2253 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2254 { 2255 PetscErrorCode ierr; 2256 PetscInt ii; 2257 Mat_SeqSBAIJ *sbaij; 2258 2259 PetscFunctionBegin; 2260 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2261 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2262 2263 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2264 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2265 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2266 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2267 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2268 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2269 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2270 2271 sbaij->i = i; 2272 sbaij->j = j; 2273 sbaij->a = a; 2274 sbaij->singlemalloc = PETSC_FALSE; 2275 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2276 sbaij->free_a = PETSC_FALSE; 2277 sbaij->free_ij = PETSC_FALSE; 2278 2279 for (ii=0; ii<m; ii++) { 2280 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2281 #if defined(PETSC_USE_DEBUG) 2282 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]); 2283 #endif 2284 } 2285 #if defined(PETSC_USE_DEBUG) 2286 for (ii=0; ii<sbaij->i[m]; ii++) { 2287 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2288 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2289 } 2290 #endif 2291 2292 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2293 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2294 PetscFunctionReturn(0); 2295 } 2296 2297 2298 2299 2300 2301