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(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,"Matrix Object");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(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(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(inA,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(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 }; 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1521 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1522 { 1523 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1524 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1525 PetscErrorCode ierr; 1526 1527 PetscFunctionBegin; 1528 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1529 1530 /* allocate space for values if not already there */ 1531 if (!aij->saved_values) { 1532 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1533 } 1534 1535 /* copy values over */ 1536 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1542 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1543 { 1544 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1545 PetscErrorCode ierr; 1546 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 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 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1551 1552 /* copy values over */ 1553 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1559 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1560 { 1561 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1562 PetscErrorCode ierr; 1563 PetscInt i,mbs,bs2; 1564 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1565 1566 PetscFunctionBegin; 1567 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1568 B->preallocated = PETSC_TRUE; 1569 1570 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1571 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1572 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1573 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1574 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1575 1576 mbs = B->rmap->N/bs; 1577 bs2 = bs*bs; 1578 1579 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1580 1581 if (nz == MAT_SKIP_ALLOCATION) { 1582 skipallocation = PETSC_TRUE; 1583 nz = 0; 1584 } 1585 1586 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1587 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1588 if (nnz) { 1589 for (i=0; i<mbs; i++) { 1590 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]); 1591 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); 1592 } 1593 } 1594 1595 B->ops->mult = MatMult_SeqSBAIJ_N; 1596 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1597 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1598 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1599 1600 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1601 if (!flg) { 1602 switch (bs) { 1603 case 1: 1604 B->ops->mult = MatMult_SeqSBAIJ_1; 1605 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1606 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1607 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1608 break; 1609 case 2: 1610 B->ops->mult = MatMult_SeqSBAIJ_2; 1611 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1612 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1613 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1614 break; 1615 case 3: 1616 B->ops->mult = MatMult_SeqSBAIJ_3; 1617 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1618 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1619 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1620 break; 1621 case 4: 1622 B->ops->mult = MatMult_SeqSBAIJ_4; 1623 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1624 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1625 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1626 break; 1627 case 5: 1628 B->ops->mult = MatMult_SeqSBAIJ_5; 1629 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1630 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1631 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1632 break; 1633 case 6: 1634 B->ops->mult = MatMult_SeqSBAIJ_6; 1635 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1636 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1637 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1638 break; 1639 case 7: 1640 B->ops->mult = MatMult_SeqSBAIJ_7; 1641 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1642 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1643 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1644 break; 1645 } 1646 } 1647 1648 b->mbs = mbs; 1649 b->nbs = mbs; 1650 if (!skipallocation) { 1651 if (!b->imax) { 1652 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1653 1654 b->free_imax_ilen = PETSC_TRUE; 1655 1656 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1657 } 1658 if (!nnz) { 1659 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1660 else if (nz <= 0) nz = 1; 1661 for (i=0; i<mbs; i++) b->imax[i] = nz; 1662 nz = nz*mbs; /* total nz */ 1663 } else { 1664 nz = 0; 1665 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1666 } 1667 /* b->ilen will count nonzeros in each block row so far. */ 1668 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1669 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1670 1671 /* allocate the matrix space */ 1672 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1673 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1674 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1675 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1676 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1677 1678 b->singlemalloc = PETSC_TRUE; 1679 1680 /* pointer to beginning of each row */ 1681 b->i[0] = 0; 1682 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1683 1684 b->free_a = PETSC_TRUE; 1685 b->free_ij = PETSC_TRUE; 1686 } else { 1687 b->free_a = PETSC_FALSE; 1688 b->free_ij = PETSC_FALSE; 1689 } 1690 1691 B->rmap->bs = bs; 1692 b->bs2 = bs2; 1693 b->nz = 0; 1694 b->maxnz = nz; 1695 1696 b->inew = 0; 1697 b->jnew = 0; 1698 b->anew = 0; 1699 b->a2anew = 0; 1700 b->permute = PETSC_FALSE; 1701 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1702 PetscFunctionReturn(0); 1703 } 1704 1705 /* 1706 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1707 */ 1708 #undef __FUNCT__ 1709 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1710 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1711 { 1712 PetscErrorCode ierr; 1713 PetscBool flg = PETSC_FALSE; 1714 PetscInt bs = B->rmap->bs; 1715 1716 PetscFunctionBegin; 1717 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1718 if (flg) bs = 8; 1719 1720 if (!natural) { 1721 switch (bs) { 1722 case 1: 1723 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1724 break; 1725 case 2: 1726 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1727 break; 1728 case 3: 1729 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1730 break; 1731 case 4: 1732 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1733 break; 1734 case 5: 1735 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1736 break; 1737 case 6: 1738 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1739 break; 1740 case 7: 1741 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1742 break; 1743 default: 1744 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1745 break; 1746 } 1747 } else { 1748 switch (bs) { 1749 case 1: 1750 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1751 break; 1752 case 2: 1753 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1754 break; 1755 case 3: 1756 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1757 break; 1758 case 4: 1759 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1760 break; 1761 case 5: 1762 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1763 break; 1764 case 6: 1765 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1766 break; 1767 case 7: 1768 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1769 break; 1770 default: 1771 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1772 break; 1773 } 1774 } 1775 PetscFunctionReturn(0); 1776 } 1777 1778 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1779 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1780 1781 #undef __FUNCT__ 1782 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1783 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1784 { 1785 PetscInt n = A->rmap->n; 1786 PetscErrorCode ierr; 1787 1788 PetscFunctionBegin; 1789 #if defined(PETSC_USE_COMPLEX) 1790 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1791 #endif 1792 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1793 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1794 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1795 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1796 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1797 1798 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1799 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1800 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1801 (*B)->factortype = ftype; 1802 PetscFunctionReturn(0); 1803 } 1804 1805 #undef __FUNCT__ 1806 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1807 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 1808 { 1809 PetscFunctionBegin; 1810 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1811 *flg = PETSC_TRUE; 1812 } else { 1813 *flg = PETSC_FALSE; 1814 } 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #if defined(PETSC_HAVE_MUMPS) 1819 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1820 #endif 1821 #if defined(PETSC_HAVE_PASTIX) 1822 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1823 #endif 1824 #if defined(PETSC_HAVE_CHOLMOD) 1825 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1826 #endif 1827 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_sbstrm(Mat,MatFactorType,Mat*); 1828 1829 /*MC 1830 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1831 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1832 1833 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1834 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1835 1836 Options Database Keys: 1837 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1838 1839 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1840 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1841 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. 1842 1843 1844 Level: beginner 1845 1846 .seealso: MatCreateSeqSBAIJ 1847 M*/ 1848 1849 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1850 1851 #undef __FUNCT__ 1852 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1853 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1854 { 1855 Mat_SeqSBAIJ *b; 1856 PetscErrorCode ierr; 1857 PetscMPIInt size; 1858 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1859 1860 PetscFunctionBegin; 1861 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1862 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1863 1864 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1865 B->data = (void*)b; 1866 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1867 1868 B->ops->destroy = MatDestroy_SeqSBAIJ; 1869 B->ops->view = MatView_SeqSBAIJ; 1870 b->row = 0; 1871 b->icol = 0; 1872 b->reallocs = 0; 1873 b->saved_values = 0; 1874 b->inode.limit = 5; 1875 b->inode.max_limit = 5; 1876 1877 b->roworiented = PETSC_TRUE; 1878 b->nonew = 0; 1879 b->diag = 0; 1880 b->solve_work = 0; 1881 b->mult_work = 0; 1882 B->spptr = 0; 1883 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1884 b->keepnonzeropattern = PETSC_FALSE; 1885 b->xtoy = 0; 1886 b->XtoY = 0; 1887 1888 b->inew = 0; 1889 b->jnew = 0; 1890 b->anew = 0; 1891 b->a2anew = 0; 1892 b->permute = PETSC_FALSE; 1893 1894 b->ignore_ltriangular = PETSC_TRUE; 1895 1896 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1897 1898 b->getrow_utriangular = PETSC_FALSE; 1899 1900 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1901 1902 #if defined(PETSC_HAVE_PASTIX) 1903 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1904 #endif 1905 #if defined(PETSC_HAVE_MUMPS) 1906 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1907 #endif 1908 #if defined(PETSC_HAVE_CHOLMOD) 1909 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1910 #endif 1911 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1912 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_sbstrm_C",MatGetFactor_seqsbaij_sbstrm);CHKERRQ(ierr); 1914 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1915 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1916 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1917 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1918 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1919 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1920 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1921 1922 B->symmetric = PETSC_TRUE; 1923 B->structurally_symmetric = PETSC_TRUE; 1924 B->symmetric_set = PETSC_TRUE; 1925 B->structurally_symmetric_set = PETSC_TRUE; 1926 1927 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1928 1929 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1930 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1931 if (no_unroll) { 1932 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1933 } 1934 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1935 if (no_inode) { 1936 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1937 } 1938 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1939 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1940 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1941 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1942 PetscFunctionReturn(0); 1943 } 1944 1945 #undef __FUNCT__ 1946 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1947 /*@C 1948 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1949 compressed row) format. For good matrix assembly performance the 1950 user should preallocate the matrix storage by setting the parameter nz 1951 (or the array nnz). By setting these parameters accurately, performance 1952 during matrix assembly can be increased by more than a factor of 50. 1953 1954 Collective on Mat 1955 1956 Input Parameters: 1957 + A - the symmetric matrix 1958 . bs - size of block 1959 . nz - number of block nonzeros per block row (same for all rows) 1960 - nnz - array containing the number of block nonzeros in the upper triangular plus 1961 diagonal portion of each block (possibly different for each block row) or NULL 1962 1963 Options Database Keys: 1964 . -mat_no_unroll - uses code that does not unroll the loops in the 1965 block calculations (much slower) 1966 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1967 1968 Level: intermediate 1969 1970 Notes: 1971 Specify the preallocated storage with either nz or nnz (not both). 1972 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 1973 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 1974 1975 You can call MatGetInfo() to get information on how effective the preallocation was; 1976 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1977 You can also run with the option -info and look for messages with the string 1978 malloc in them to see if additional memory allocation was needed. 1979 1980 If the nnz parameter is given then the nz parameter is ignored 1981 1982 1983 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 1984 @*/ 1985 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1986 { 1987 PetscErrorCode ierr; 1988 1989 PetscFunctionBegin; 1990 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1991 PetscValidType(B,1); 1992 PetscValidLogicalCollectiveInt(B,bs,2); 1993 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 1994 PetscFunctionReturn(0); 1995 } 1996 1997 #undef __FUNCT__ 1998 #define __FUNCT__ "MatCreateSeqSBAIJ" 1999 /*@C 2000 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2001 compressed row) format. For good matrix assembly performance the 2002 user should preallocate the matrix storage by setting the parameter nz 2003 (or the array nnz). By setting these parameters accurately, performance 2004 during matrix assembly can be increased by more than a factor of 50. 2005 2006 Collective on MPI_Comm 2007 2008 Input Parameters: 2009 + comm - MPI communicator, set to PETSC_COMM_SELF 2010 . bs - size of block 2011 . m - number of rows, or number of columns 2012 . nz - number of block nonzeros per block row (same for all rows) 2013 - nnz - array containing the number of block nonzeros in the upper triangular plus 2014 diagonal portion of each block (possibly different for each block row) or NULL 2015 2016 Output Parameter: 2017 . A - the symmetric matrix 2018 2019 Options Database Keys: 2020 . -mat_no_unroll - uses code that does not unroll the loops in the 2021 block calculations (much slower) 2022 . -mat_block_size - size of the blocks to use 2023 2024 Level: intermediate 2025 2026 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2027 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2028 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2029 2030 Notes: 2031 The number of rows and columns must be divisible by blocksize. 2032 This matrix type does not support complex Hermitian operation. 2033 2034 Specify the preallocated storage with either nz or nnz (not both). 2035 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2036 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2037 2038 If the nnz parameter is given then the nz parameter is ignored 2039 2040 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2041 @*/ 2042 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2043 { 2044 PetscErrorCode ierr; 2045 2046 PetscFunctionBegin; 2047 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2048 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2049 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2050 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 #undef __FUNCT__ 2055 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2056 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2057 { 2058 Mat C; 2059 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2060 PetscErrorCode ierr; 2061 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2062 2063 PetscFunctionBegin; 2064 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2065 2066 *B = 0; 2067 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2068 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2069 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2070 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2071 c = (Mat_SeqSBAIJ*)C->data; 2072 2073 C->preallocated = PETSC_TRUE; 2074 C->factortype = A->factortype; 2075 c->row = 0; 2076 c->icol = 0; 2077 c->saved_values = 0; 2078 c->keepnonzeropattern = a->keepnonzeropattern; 2079 C->assembled = PETSC_TRUE; 2080 2081 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2082 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2083 c->bs2 = a->bs2; 2084 c->mbs = a->mbs; 2085 c->nbs = a->nbs; 2086 2087 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2088 c->imax = a->imax; 2089 c->ilen = a->ilen; 2090 c->free_imax_ilen = PETSC_FALSE; 2091 } else { 2092 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2093 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2094 for (i=0; i<mbs; i++) { 2095 c->imax[i] = a->imax[i]; 2096 c->ilen[i] = a->ilen[i]; 2097 } 2098 c->free_imax_ilen = PETSC_TRUE; 2099 } 2100 2101 /* allocate the matrix space */ 2102 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2103 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2104 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2105 c->i = a->i; 2106 c->j = a->j; 2107 c->singlemalloc = PETSC_FALSE; 2108 c->free_a = PETSC_TRUE; 2109 c->free_ij = PETSC_FALSE; 2110 c->parent = A; 2111 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2112 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2113 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2114 } else { 2115 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2116 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2117 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2118 c->singlemalloc = PETSC_TRUE; 2119 c->free_a = PETSC_TRUE; 2120 c->free_ij = PETSC_TRUE; 2121 } 2122 if (mbs > 0) { 2123 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2124 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2125 } 2126 if (cpvalues == MAT_COPY_VALUES) { 2127 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2128 } else { 2129 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2130 } 2131 if (a->jshort) { 2132 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2133 /* if the parent matrix is reassembled, this child matrix will never notice */ 2134 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2135 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2136 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2137 2138 c->free_jshort = PETSC_TRUE; 2139 } 2140 } 2141 2142 c->roworiented = a->roworiented; 2143 c->nonew = a->nonew; 2144 2145 if (a->diag) { 2146 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2147 c->diag = a->diag; 2148 c->free_diag = PETSC_FALSE; 2149 } else { 2150 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2151 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2152 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2153 c->free_diag = PETSC_TRUE; 2154 } 2155 } 2156 c->nz = a->nz; 2157 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2158 c->solve_work = 0; 2159 c->mult_work = 0; 2160 2161 *B = C; 2162 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2163 PetscFunctionReturn(0); 2164 } 2165 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2168 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2169 { 2170 Mat_SeqSBAIJ *a; 2171 PetscErrorCode ierr; 2172 int fd; 2173 PetscMPIInt size; 2174 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2175 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2176 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2177 PetscInt *masked,nmask,tmp,bs2,ishift; 2178 PetscScalar *aa; 2179 MPI_Comm comm; 2180 2181 PetscFunctionBegin; 2182 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2183 ierr = PetscOptionsGetInt(((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2184 bs2 = bs*bs; 2185 2186 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2187 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2188 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2189 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2190 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2191 M = header[1]; N = header[2]; nz = header[3]; 2192 2193 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2194 2195 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2196 2197 /* 2198 This code adds extra rows to make sure the number of rows is 2199 divisible by the blocksize 2200 */ 2201 mbs = M/bs; 2202 extra_rows = bs - M + bs*(mbs); 2203 if (extra_rows == bs) extra_rows = 0; 2204 else mbs++; 2205 if (extra_rows) { 2206 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2207 } 2208 2209 /* Set global sizes if not already set */ 2210 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2211 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2212 } else { /* Check if the matrix global sizes are correct */ 2213 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2214 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); 2215 } 2216 2217 /* read in row lengths */ 2218 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2219 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2220 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2221 2222 /* read in column indices */ 2223 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2224 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2225 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2226 2227 /* loop over row lengths determining block row lengths */ 2228 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2229 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2230 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2231 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2232 rowcount = 0; 2233 nzcount = 0; 2234 for (i=0; i<mbs; i++) { 2235 nmask = 0; 2236 for (j=0; j<bs; j++) { 2237 kmax = rowlengths[rowcount]; 2238 for (k=0; k<kmax; k++) { 2239 tmp = jj[nzcount++]/bs; /* block col. index */ 2240 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2241 } 2242 rowcount++; 2243 } 2244 s_browlengths[i] += nmask; 2245 2246 /* zero out the mask elements we set */ 2247 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2248 } 2249 2250 /* Do preallocation */ 2251 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2252 a = (Mat_SeqSBAIJ*)newmat->data; 2253 2254 /* set matrix "i" values */ 2255 a->i[0] = 0; 2256 for (i=1; i<= mbs; i++) { 2257 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2258 a->ilen[i-1] = s_browlengths[i-1]; 2259 } 2260 a->nz = a->i[mbs]; 2261 2262 /* read in nonzero values */ 2263 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2264 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2265 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2266 2267 /* set "a" and "j" values into matrix */ 2268 nzcount = 0; jcount = 0; 2269 for (i=0; i<mbs; i++) { 2270 nzcountb = nzcount; 2271 nmask = 0; 2272 for (j=0; j<bs; j++) { 2273 kmax = rowlengths[i*bs+j]; 2274 for (k=0; k<kmax; k++) { 2275 tmp = jj[nzcount++]/bs; /* block col. index */ 2276 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2277 } 2278 } 2279 /* sort the masked values */ 2280 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2281 2282 /* set "j" values into matrix */ 2283 maskcount = 1; 2284 for (j=0; j<nmask; j++) { 2285 a->j[jcount++] = masked[j]; 2286 mask[masked[j]] = maskcount++; 2287 } 2288 2289 /* set "a" values into matrix */ 2290 ishift = bs2*a->i[i]; 2291 for (j=0; j<bs; j++) { 2292 kmax = rowlengths[i*bs+j]; 2293 for (k=0; k<kmax; k++) { 2294 tmp = jj[nzcountb]/bs; /* block col. index */ 2295 if (tmp >= i) { 2296 block = mask[tmp] - 1; 2297 point = jj[nzcountb] - bs*tmp; 2298 idx = ishift + bs2*block + j + bs*point; 2299 a->a[idx] = aa[nzcountb]; 2300 } 2301 nzcountb++; 2302 } 2303 } 2304 /* zero out the mask elements we set */ 2305 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2306 } 2307 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2308 2309 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2310 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2311 ierr = PetscFree(aa);CHKERRQ(ierr); 2312 ierr = PetscFree(jj);CHKERRQ(ierr); 2313 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2314 2315 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2316 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 #undef __FUNCT__ 2321 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2322 /*@ 2323 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2324 (upper triangular entries in CSR format) provided by the user. 2325 2326 Collective on MPI_Comm 2327 2328 Input Parameters: 2329 + comm - must be an MPI communicator of size 1 2330 . bs - size of block 2331 . m - number of rows 2332 . n - number of columns 2333 . i - row indices 2334 . j - column indices 2335 - a - matrix values 2336 2337 Output Parameter: 2338 . mat - the matrix 2339 2340 Level: advanced 2341 2342 Notes: 2343 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2344 once the matrix is destroyed 2345 2346 You cannot set new nonzero locations into this matrix, that will generate an error. 2347 2348 The i and j indices are 0 based 2349 2350 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 2351 it is the regular CSR format excluding the lower triangular elements. 2352 2353 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2354 2355 @*/ 2356 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2357 { 2358 PetscErrorCode ierr; 2359 PetscInt ii; 2360 Mat_SeqSBAIJ *sbaij; 2361 2362 PetscFunctionBegin; 2363 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2364 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2365 2366 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2367 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2368 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2369 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2370 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2371 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2372 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2373 2374 sbaij->i = i; 2375 sbaij->j = j; 2376 sbaij->a = a; 2377 2378 sbaij->singlemalloc = PETSC_FALSE; 2379 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2380 sbaij->free_a = PETSC_FALSE; 2381 sbaij->free_ij = PETSC_FALSE; 2382 2383 for (ii=0; ii<m; ii++) { 2384 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2385 #if defined(PETSC_USE_DEBUG) 2386 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]); 2387 #endif 2388 } 2389 #if defined(PETSC_USE_DEBUG) 2390 for (ii=0; ii<sbaij->i[m]; ii++) { 2391 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2392 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]); 2393 } 2394 #endif 2395 2396 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2397 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2398 PetscFunctionReturn(0); 2399 } 2400 2401 2402 2403 2404 2405