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