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