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