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