1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 8 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 9 #include <petscblaslapack.h> 10 #include <petscbt.h> 11 #include <petsc/private/kernels/blocktranspose.h> 12 13 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 14 { 15 PetscErrorCode ierr; 16 PetscBool flg; 17 char type[256]; 18 19 PetscFunctionBegin; 20 ierr = PetscObjectOptionsBegin((PetscObject)A); 21 ierr = PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg);CHKERRQ(ierr); 22 if (flg) { 23 ierr = MatSeqAIJSetType(A,type);CHKERRQ(ierr); 24 } 25 ierr = PetscOptionsEnd();CHKERRQ(ierr); 26 PetscFunctionReturn(0); 27 } 28 29 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 30 { 31 PetscErrorCode ierr; 32 PetscInt i,m,n; 33 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 34 35 PetscFunctionBegin; 36 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 37 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 38 if (type == NORM_2) { 39 for (i=0; i<aij->i[m]; i++) { 40 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 41 } 42 } else if (type == NORM_1) { 43 for (i=0; i<aij->i[m]; i++) { 44 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 45 } 46 } else if (type == NORM_INFINITY) { 47 for (i=0; i<aij->i[m]; i++) { 48 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 49 } 50 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 51 52 if (type == NORM_2) { 53 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 54 } 55 PetscFunctionReturn(0); 56 } 57 58 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 59 { 60 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 61 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 62 const PetscInt *jj = a->j,*ii = a->i; 63 PetscInt *rows; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 for (i=0; i<m; i++) { 68 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 69 cnt++; 70 } 71 } 72 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 73 cnt = 0; 74 for (i=0; i<m; i++) { 75 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 76 rows[cnt] = i; 77 cnt++; 78 } 79 } 80 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 85 { 86 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 87 const MatScalar *aa = a->a; 88 PetscInt i,m=A->rmap->n,cnt = 0; 89 const PetscInt *ii = a->i,*jj = a->j,*diag; 90 PetscInt *rows; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 95 diag = a->diag; 96 for (i=0; i<m; i++) { 97 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 98 cnt++; 99 } 100 } 101 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 102 cnt = 0; 103 for (i=0; i<m; i++) { 104 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 105 rows[cnt++] = i; 106 } 107 } 108 *nrows = cnt; 109 *zrows = rows; 110 PetscFunctionReturn(0); 111 } 112 113 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 114 { 115 PetscInt nrows,*rows; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 *zrows = NULL; 120 ierr = MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);CHKERRQ(ierr); 121 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 126 { 127 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 128 const MatScalar *aa; 129 PetscInt m=A->rmap->n,cnt = 0; 130 const PetscInt *ii; 131 PetscInt n,i,j,*rows; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 *keptrows = 0; 136 ii = a->i; 137 for (i=0; i<m; i++) { 138 n = ii[i+1] - ii[i]; 139 if (!n) { 140 cnt++; 141 goto ok1; 142 } 143 aa = a->a + ii[i]; 144 for (j=0; j<n; j++) { 145 if (aa[j] != 0.0) goto ok1; 146 } 147 cnt++; 148 ok1:; 149 } 150 if (!cnt) PetscFunctionReturn(0); 151 ierr = PetscMalloc1(A->rmap->n-cnt,&rows);CHKERRQ(ierr); 152 cnt = 0; 153 for (i=0; i<m; i++) { 154 n = ii[i+1] - ii[i]; 155 if (!n) continue; 156 aa = a->a + ii[i]; 157 for (j=0; j<n; j++) { 158 if (aa[j] != 0.0) { 159 rows[cnt++] = i; 160 break; 161 } 162 } 163 } 164 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 169 { 170 PetscErrorCode ierr; 171 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 172 PetscInt i,m = Y->rmap->n; 173 const PetscInt *diag; 174 MatScalar *aa = aij->a; 175 const PetscScalar *v; 176 PetscBool missing; 177 178 PetscFunctionBegin; 179 if (Y->assembled) { 180 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);CHKERRQ(ierr); 181 if (!missing) { 182 diag = aij->diag; 183 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 184 if (is == INSERT_VALUES) { 185 for (i=0; i<m; i++) { 186 aa[diag[i]] = v[i]; 187 } 188 } else { 189 for (i=0; i<m; i++) { 190 aa[diag[i]] += v[i]; 191 } 192 } 193 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 197 } 198 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 203 { 204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 205 PetscErrorCode ierr; 206 PetscInt i,ishift; 207 208 PetscFunctionBegin; 209 *m = A->rmap->n; 210 if (!ia) PetscFunctionReturn(0); 211 ishift = 0; 212 if (symmetric && !A->structurally_symmetric) { 213 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 214 } else if (oshift == 1) { 215 PetscInt *tia; 216 PetscInt nz = a->i[A->rmap->n]; 217 /* malloc space and add 1 to i and j indices */ 218 ierr = PetscMalloc1(A->rmap->n+1,&tia);CHKERRQ(ierr); 219 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 220 *ia = tia; 221 if (ja) { 222 PetscInt *tja; 223 ierr = PetscMalloc1(nz+1,&tja);CHKERRQ(ierr); 224 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 225 *ja = tja; 226 } 227 } else { 228 *ia = a->i; 229 if (ja) *ja = a->j; 230 } 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 if (!ia) PetscFunctionReturn(0); 240 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 241 ierr = PetscFree(*ia);CHKERRQ(ierr); 242 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 243 } 244 PetscFunctionReturn(0); 245 } 246 247 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 248 { 249 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 250 PetscErrorCode ierr; 251 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 252 PetscInt nz = a->i[m],row,*jj,mr,col; 253 254 PetscFunctionBegin; 255 *nn = n; 256 if (!ia) PetscFunctionReturn(0); 257 if (symmetric) { 258 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 259 } else { 260 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 261 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 262 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 263 jj = a->j; 264 for (i=0; i<nz; i++) { 265 collengths[jj[i]]++; 266 } 267 cia[0] = oshift; 268 for (i=0; i<n; i++) { 269 cia[i+1] = cia[i] + collengths[i]; 270 } 271 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 272 jj = a->j; 273 for (row=0; row<m; row++) { 274 mr = a->i[row+1] - a->i[row]; 275 for (i=0; i<mr; i++) { 276 col = *jj++; 277 278 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 279 } 280 } 281 ierr = PetscFree(collengths);CHKERRQ(ierr); 282 *ia = cia; *ja = cja; 283 } 284 PetscFunctionReturn(0); 285 } 286 287 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 288 { 289 PetscErrorCode ierr; 290 291 PetscFunctionBegin; 292 if (!ia) PetscFunctionReturn(0); 293 294 ierr = PetscFree(*ia);CHKERRQ(ierr); 295 ierr = PetscFree(*ja);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 301 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 302 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 303 */ 304 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 305 { 306 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 307 PetscErrorCode ierr; 308 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 309 PetscInt nz = a->i[m],row,*jj,mr,col; 310 PetscInt *cspidx; 311 312 PetscFunctionBegin; 313 *nn = n; 314 if (!ia) PetscFunctionReturn(0); 315 316 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 317 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 318 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 319 ierr = PetscMalloc1(nz+1,&cspidx);CHKERRQ(ierr); 320 jj = a->j; 321 for (i=0; i<nz; i++) { 322 collengths[jj[i]]++; 323 } 324 cia[0] = oshift; 325 for (i=0; i<n; i++) { 326 cia[i+1] = cia[i] + collengths[i]; 327 } 328 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 329 jj = a->j; 330 for (row=0; row<m; row++) { 331 mr = a->i[row+1] - a->i[row]; 332 for (i=0; i<mr; i++) { 333 col = *jj++; 334 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 335 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 336 } 337 } 338 ierr = PetscFree(collengths);CHKERRQ(ierr); 339 *ia = cia; *ja = cja; 340 *spidx = cspidx; 341 PetscFunctionReturn(0); 342 } 343 344 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 345 { 346 PetscErrorCode ierr; 347 348 PetscFunctionBegin; 349 ierr = MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);CHKERRQ(ierr); 350 ierr = PetscFree(*spidx);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 355 { 356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 357 PetscInt *ai = a->i; 358 PetscErrorCode ierr; 359 360 PetscFunctionBegin; 361 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 362 PetscFunctionReturn(0); 363 } 364 365 /* 366 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 367 368 - a single row of values is set with each call 369 - no row or column indices are negative or (in error) larger than the number of rows or columns 370 - the values are always added to the matrix, not set 371 - no new locations are introduced in the nonzero structure of the matrix 372 373 This does NOT assume the global column indices are sorted 374 375 */ 376 377 #include <petsc/private/isimpl.h> 378 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 379 { 380 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 381 PetscInt low,high,t,row,nrow,i,col,l; 382 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 383 PetscInt lastcol = -1; 384 MatScalar *ap,value,*aa = a->a; 385 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 386 387 row = ridx[im[0]]; 388 rp = aj + ai[row]; 389 ap = aa + ai[row]; 390 nrow = ailen[row]; 391 low = 0; 392 high = nrow; 393 for (l=0; l<n; l++) { /* loop over added columns */ 394 col = cidx[in[l]]; 395 value = v[l]; 396 397 if (col <= lastcol) low = 0; 398 else high = nrow; 399 lastcol = col; 400 while (high-low > 5) { 401 t = (low+high)/2; 402 if (rp[t] > col) high = t; 403 else low = t; 404 } 405 for (i=low; i<high; i++) { 406 if (rp[i] == col) { 407 ap[i] += value; 408 low = i + 1; 409 break; 410 } 411 } 412 } 413 return 0; 414 } 415 416 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 417 { 418 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 419 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 420 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 421 PetscErrorCode ierr; 422 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 423 MatScalar *ap=NULL,value=0.0,*aa = a->a; 424 PetscBool ignorezeroentries = a->ignorezeroentries; 425 PetscBool roworiented = a->roworiented; 426 427 PetscFunctionBegin; 428 for (k=0; k<m; k++) { /* loop over added rows */ 429 row = im[k]; 430 if (row < 0) continue; 431 #if defined(PETSC_USE_DEBUG) 432 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); 433 #endif 434 rp = aj + ai[row]; 435 if (!A->structure_only) ap = aa + ai[row]; 436 rmax = imax[row]; nrow = ailen[row]; 437 low = 0; 438 high = nrow; 439 for (l=0; l<n; l++) { /* loop over added columns */ 440 if (in[l] < 0) continue; 441 #if defined(PETSC_USE_DEBUG) 442 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); 443 #endif 444 col = in[l]; 445 if (!A->structure_only) { 446 if (roworiented) { 447 value = v[l + k*n]; 448 } else { 449 value = v[k + l*m]; 450 } 451 } else { /* A->structure_only */ 452 value = 1; /* avoid 'continue' below? */ 453 } 454 if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES) && row != col) continue; 455 456 if (col <= lastcol) low = 0; 457 else high = nrow; 458 lastcol = col; 459 while (high-low > 5) { 460 t = (low+high)/2; 461 if (rp[t] > col) high = t; 462 else low = t; 463 } 464 for (i=low; i<high; i++) { 465 if (rp[i] > col) break; 466 if (rp[i] == col) { 467 if (!A->structure_only) { 468 if (is == ADD_VALUES) ap[i] += value; 469 else ap[i] = value; 470 } 471 low = i + 1; 472 goto noinsert; 473 } 474 } 475 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 476 if (nonew == 1) goto noinsert; 477 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 478 if (A->structure_only) { 479 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 480 } else { 481 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 482 } 483 N = nrow++ - 1; a->nz++; high++; 484 /* shift up all the later entries in this row */ 485 for (ii=N; ii>=i; ii--) { 486 rp[ii+1] = rp[ii]; 487 if (!A->structure_only) ap[ii+1] = ap[ii]; 488 } 489 rp[i] = col; 490 if (!A->structure_only) ap[i] = value; 491 low = i + 1; 492 A->nonzerostate++; 493 noinsert:; 494 } 495 ailen[row] = nrow; 496 } 497 PetscFunctionReturn(0); 498 } 499 500 501 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 502 { 503 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 504 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 505 PetscInt *ai = a->i,*ailen = a->ilen; 506 MatScalar *ap,*aa = a->a; 507 508 PetscFunctionBegin; 509 for (k=0; k<m; k++) { /* loop over rows */ 510 row = im[k]; 511 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 512 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); 513 rp = aj + ai[row]; ap = aa + ai[row]; 514 nrow = ailen[row]; 515 for (l=0; l<n; l++) { /* loop over columns */ 516 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 517 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); 518 col = in[l]; 519 high = nrow; low = 0; /* assume unsorted */ 520 while (high-low > 5) { 521 t = (low+high)/2; 522 if (rp[t] > col) high = t; 523 else low = t; 524 } 525 for (i=low; i<high; i++) { 526 if (rp[i] > col) break; 527 if (rp[i] == col) { 528 *v++ = ap[i]; 529 goto finished; 530 } 531 } 532 *v++ = 0.0; 533 finished:; 534 } 535 } 536 PetscFunctionReturn(0); 537 } 538 539 540 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 541 { 542 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 543 PetscErrorCode ierr; 544 PetscInt i,*col_lens; 545 int fd; 546 FILE *file; 547 548 PetscFunctionBegin; 549 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 550 ierr = PetscMalloc1(4+A->rmap->n,&col_lens);CHKERRQ(ierr); 551 552 col_lens[0] = MAT_FILE_CLASSID; 553 col_lens[1] = A->rmap->n; 554 col_lens[2] = A->cmap->n; 555 col_lens[3] = a->nz; 556 557 /* store lengths of each row and write (including header) to file */ 558 for (i=0; i<A->rmap->n; i++) { 559 col_lens[4+i] = a->i[i+1] - a->i[i]; 560 } 561 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 562 ierr = PetscFree(col_lens);CHKERRQ(ierr); 563 564 /* store column indices (zero start index) */ 565 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 566 567 /* store nonzero values */ 568 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 569 570 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 571 if (file) { 572 fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs)); 573 } 574 PetscFunctionReturn(0); 575 } 576 577 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 578 { 579 PetscErrorCode ierr; 580 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 581 PetscInt i,k,m=A->rmap->N; 582 583 PetscFunctionBegin; 584 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 585 for (i=0; i<m; i++) { 586 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 587 for (k=a->i[i]; k<a->i[i+1]; k++) { 588 ierr = PetscViewerASCIIPrintf(viewer," (%D) ",a->j[k]);CHKERRQ(ierr); 589 } 590 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 591 } 592 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 597 598 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 599 { 600 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 601 PetscErrorCode ierr; 602 PetscInt i,j,m = A->rmap->n; 603 const char *name; 604 PetscViewerFormat format; 605 606 PetscFunctionBegin; 607 if (A->structure_only) { 608 ierr = MatView_SeqAIJ_ASCII_structonly(A,viewer);CHKERRQ(ierr); 609 PetscFunctionReturn(0); 610 } 611 612 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 613 if (format == PETSC_VIEWER_ASCII_MATLAB) { 614 PetscInt nofinalvalue = 0; 615 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 616 /* Need a dummy value to ensure the dimension of the matrix. */ 617 nofinalvalue = 1; 618 } 619 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 620 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 621 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 622 #if defined(PETSC_USE_COMPLEX) 623 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,4);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 624 #else 625 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 626 #endif 627 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 628 629 for (i=0; i<m; i++) { 630 for (j=a->i[i]; j<a->i[i+1]; j++) { 631 #if defined(PETSC_USE_COMPLEX) 632 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 633 #else 634 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);CHKERRQ(ierr); 635 #endif 636 } 637 } 638 if (nofinalvalue) { 639 #if defined(PETSC_USE_COMPLEX) 640 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",m,A->cmap->n,0.,0.);CHKERRQ(ierr); 641 #else 642 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 643 #endif 644 } 645 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 646 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 647 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 648 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 649 PetscFunctionReturn(0); 650 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 651 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 652 for (i=0; i<m; i++) { 653 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 654 for (j=a->i[i]; j<a->i[i+1]; j++) { 655 #if defined(PETSC_USE_COMPLEX) 656 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 657 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 658 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 659 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 660 } else if (PetscRealPart(a->a[j]) != 0.0) { 661 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 662 } 663 #else 664 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr);} 665 #endif 666 } 667 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 668 } 669 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 670 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 671 PetscInt nzd=0,fshift=1,*sptr; 672 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 673 ierr = PetscMalloc1(m+1,&sptr);CHKERRQ(ierr); 674 for (i=0; i<m; i++) { 675 sptr[i] = nzd+1; 676 for (j=a->i[i]; j<a->i[i+1]; j++) { 677 if (a->j[j] >= i) { 678 #if defined(PETSC_USE_COMPLEX) 679 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 680 #else 681 if (a->a[j] != 0.0) nzd++; 682 #endif 683 } 684 } 685 } 686 sptr[m] = nzd+1; 687 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 688 for (i=0; i<m+1; i+=6) { 689 if (i+4<m) { 690 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr); 691 } else if (i+3<m) { 692 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr); 693 } else if (i+2<m) { 694 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr); 695 } else if (i+1<m) { 696 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr); 697 } else if (i<m) { 698 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr); 699 } else { 700 ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr); 701 } 702 } 703 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 704 ierr = PetscFree(sptr);CHKERRQ(ierr); 705 for (i=0; i<m; i++) { 706 for (j=a->i[i]; j<a->i[i+1]; j++) { 707 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 708 } 709 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 710 } 711 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 712 for (i=0; i<m; i++) { 713 for (j=a->i[i]; j<a->i[i+1]; j++) { 714 if (a->j[j] >= i) { 715 #if defined(PETSC_USE_COMPLEX) 716 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 717 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 718 } 719 #else 720 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);CHKERRQ(ierr);} 721 #endif 722 } 723 } 724 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 725 } 726 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 727 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 728 PetscInt cnt = 0,jcnt; 729 PetscScalar value; 730 #if defined(PETSC_USE_COMPLEX) 731 PetscBool realonly = PETSC_TRUE; 732 733 for (i=0; i<a->i[m]; i++) { 734 if (PetscImaginaryPart(a->a[i]) != 0.0) { 735 realonly = PETSC_FALSE; 736 break; 737 } 738 } 739 #endif 740 741 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 742 for (i=0; i<m; i++) { 743 jcnt = 0; 744 for (j=0; j<A->cmap->n; j++) { 745 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 746 value = a->a[cnt++]; 747 jcnt++; 748 } else { 749 value = 0.0; 750 } 751 #if defined(PETSC_USE_COMPLEX) 752 if (realonly) { 753 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));CHKERRQ(ierr); 754 } else { 755 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));CHKERRQ(ierr); 756 } 757 #else 758 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);CHKERRQ(ierr); 759 #endif 760 } 761 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 762 } 763 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 764 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 765 PetscInt fshift=1; 766 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 767 #if defined(PETSC_USE_COMPLEX) 768 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n");CHKERRQ(ierr); 769 #else 770 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n");CHKERRQ(ierr); 771 #endif 772 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 773 for (i=0; i<m; i++) { 774 for (j=a->i[i]; j<a->i[i+1]; j++) { 775 #if defined(PETSC_USE_COMPLEX) 776 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 777 #else 778 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);CHKERRQ(ierr); 779 #endif 780 } 781 } 782 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 783 } else { 784 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 785 if (A->factortype) { 786 for (i=0; i<m; i++) { 787 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 788 /* L part */ 789 for (j=a->i[i]; j<a->i[i+1]; j++) { 790 #if defined(PETSC_USE_COMPLEX) 791 if (PetscImaginaryPart(a->a[j]) > 0.0) { 792 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 793 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 794 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 795 } else { 796 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 797 } 798 #else 799 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 800 #endif 801 } 802 /* diagonal */ 803 j = a->diag[i]; 804 #if defined(PETSC_USE_COMPLEX) 805 if (PetscImaginaryPart(a->a[j]) > 0.0) { 806 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]));CHKERRQ(ierr); 807 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 808 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])));CHKERRQ(ierr); 809 } else { 810 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 811 } 812 #else 813 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)(1.0/a->a[j]));CHKERRQ(ierr); 814 #endif 815 816 /* U part */ 817 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 818 #if defined(PETSC_USE_COMPLEX) 819 if (PetscImaginaryPart(a->a[j]) > 0.0) { 820 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 821 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 822 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 823 } else { 824 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 825 } 826 #else 827 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 828 #endif 829 } 830 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 831 } 832 } else { 833 for (i=0; i<m; i++) { 834 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 835 for (j=a->i[i]; j<a->i[i+1]; j++) { 836 #if defined(PETSC_USE_COMPLEX) 837 if (PetscImaginaryPart(a->a[j]) > 0.0) { 838 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 839 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 840 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 841 } else { 842 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 843 } 844 #else 845 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 846 #endif 847 } 848 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 849 } 850 } 851 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 852 } 853 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 854 PetscFunctionReturn(0); 855 } 856 857 #include <petscdraw.h> 858 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 859 { 860 Mat A = (Mat) Aa; 861 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 862 PetscErrorCode ierr; 863 PetscInt i,j,m = A->rmap->n; 864 int color; 865 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 866 PetscViewer viewer; 867 PetscViewerFormat format; 868 869 PetscFunctionBegin; 870 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 871 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 872 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 873 874 /* loop over matrix elements drawing boxes */ 875 876 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 877 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 878 /* Blue for negative, Cyan for zero and Red for positive */ 879 color = PETSC_DRAW_BLUE; 880 for (i=0; i<m; i++) { 881 y_l = m - i - 1.0; y_r = y_l + 1.0; 882 for (j=a->i[i]; j<a->i[i+1]; j++) { 883 x_l = a->j[j]; x_r = x_l + 1.0; 884 if (PetscRealPart(a->a[j]) >= 0.) continue; 885 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 886 } 887 } 888 color = PETSC_DRAW_CYAN; 889 for (i=0; i<m; i++) { 890 y_l = m - i - 1.0; y_r = y_l + 1.0; 891 for (j=a->i[i]; j<a->i[i+1]; j++) { 892 x_l = a->j[j]; x_r = x_l + 1.0; 893 if (a->a[j] != 0.) continue; 894 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 895 } 896 } 897 color = PETSC_DRAW_RED; 898 for (i=0; i<m; i++) { 899 y_l = m - i - 1.0; y_r = y_l + 1.0; 900 for (j=a->i[i]; j<a->i[i+1]; j++) { 901 x_l = a->j[j]; x_r = x_l + 1.0; 902 if (PetscRealPart(a->a[j]) <= 0.) continue; 903 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 904 } 905 } 906 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 907 } else { 908 /* use contour shading to indicate magnitude of values */ 909 /* first determine max of all nonzero values */ 910 PetscReal minv = 0.0, maxv = 0.0; 911 PetscInt nz = a->nz, count = 0; 912 PetscDraw popup; 913 914 for (i=0; i<nz; i++) { 915 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 916 } 917 if (minv >= maxv) maxv = minv + PETSC_SMALL; 918 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 919 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 920 921 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 922 for (i=0; i<m; i++) { 923 y_l = m - i - 1.0; 924 y_r = y_l + 1.0; 925 for (j=a->i[i]; j<a->i[i+1]; j++) { 926 x_l = a->j[j]; 927 x_r = x_l + 1.0; 928 color = PetscDrawRealToColor(PetscAbsScalar(a->a[count]),minv,maxv); 929 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 930 count++; 931 } 932 } 933 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 934 } 935 PetscFunctionReturn(0); 936 } 937 938 #include <petscdraw.h> 939 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 940 { 941 PetscErrorCode ierr; 942 PetscDraw draw; 943 PetscReal xr,yr,xl,yl,h,w; 944 PetscBool isnull; 945 946 PetscFunctionBegin; 947 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 948 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 949 if (isnull) PetscFunctionReturn(0); 950 951 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 952 xr += w; yr += h; xl = -w; yl = -h; 953 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 954 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 955 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 956 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 957 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 958 PetscFunctionReturn(0); 959 } 960 961 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 962 { 963 PetscErrorCode ierr; 964 PetscBool iascii,isbinary,isdraw; 965 966 PetscFunctionBegin; 967 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 968 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 969 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 970 if (iascii) { 971 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 972 } else if (isbinary) { 973 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 974 } else if (isdraw) { 975 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 976 } 977 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 978 PetscFunctionReturn(0); 979 } 980 981 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 982 { 983 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 984 PetscErrorCode ierr; 985 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 986 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 987 MatScalar *aa = a->a,*ap; 988 PetscReal ratio = 0.6; 989 990 PetscFunctionBegin; 991 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 992 993 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 994 for (i=1; i<m; i++) { 995 /* move each row back by the amount of empty slots (fshift) before it*/ 996 fshift += imax[i-1] - ailen[i-1]; 997 rmax = PetscMax(rmax,ailen[i]); 998 if (fshift) { 999 ip = aj + ai[i]; 1000 ap = aa + ai[i]; 1001 N = ailen[i]; 1002 for (j=0; j<N; j++) { 1003 ip[j-fshift] = ip[j]; 1004 if (!A->structure_only) ap[j-fshift] = ap[j]; 1005 } 1006 } 1007 ai[i] = ai[i-1] + ailen[i-1]; 1008 } 1009 if (m) { 1010 fshift += imax[m-1] - ailen[m-1]; 1011 ai[m] = ai[m-1] + ailen[m-1]; 1012 } 1013 1014 /* reset ilen and imax for each row */ 1015 a->nonzerorowcnt = 0; 1016 if (A->structure_only) { 1017 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1018 } else { /* !A->structure_only */ 1019 for (i=0; i<m; i++) { 1020 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1021 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1022 } 1023 } 1024 a->nz = ai[m]; 1025 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 1026 1027 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1028 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 1029 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 1030 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 1031 1032 A->info.mallocs += a->reallocs; 1033 a->reallocs = 0; 1034 A->info.nz_unneeded = (PetscReal)fshift; 1035 a->rmax = rmax; 1036 1037 if (!A->structure_only) { 1038 ierr = MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 1039 } 1040 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 1041 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1042 PetscFunctionReturn(0); 1043 } 1044 1045 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1046 { 1047 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1048 PetscInt i,nz = a->nz; 1049 MatScalar *aa = a->a; 1050 PetscErrorCode ierr; 1051 1052 PetscFunctionBegin; 1053 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1054 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1059 { 1060 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1061 PetscInt i,nz = a->nz; 1062 MatScalar *aa = a->a; 1063 PetscErrorCode ierr; 1064 1065 PetscFunctionBegin; 1066 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1067 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1072 { 1073 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1078 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1079 PetscFunctionReturn(0); 1080 } 1081 1082 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1083 { 1084 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1085 PetscErrorCode ierr; 1086 1087 PetscFunctionBegin; 1088 #if defined(PETSC_USE_LOG) 1089 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1090 #endif 1091 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1092 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1093 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1094 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1095 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1096 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1097 ierr = PetscFree(a->ipre);CHKERRQ(ierr); 1098 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1099 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1100 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1101 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1102 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1103 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1104 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1105 1106 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1107 ierr = PetscFree(A->data);CHKERRQ(ierr); 1108 1109 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1110 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1111 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1112 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1113 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1114 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1115 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1116 #if defined(PETSC_HAVE_ELEMENTAL) 1117 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL);CHKERRQ(ierr); 1118 #endif 1119 #if defined(PETSC_HAVE_HYPRE) 1120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL);CHKERRQ(ierr); 1121 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMatMult_transpose_seqaij_seqaij_C",NULL);CHKERRQ(ierr); 1122 #endif 1123 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1124 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL);CHKERRQ(ierr); 1125 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL);CHKERRQ(ierr); 1126 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1127 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1128 ierr = PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL);CHKERRQ(ierr); 1129 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1130 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1131 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_seqaij_C",NULL);CHKERRQ(ierr); 1132 PetscFunctionReturn(0); 1133 } 1134 1135 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1136 { 1137 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 switch (op) { 1142 case MAT_ROW_ORIENTED: 1143 a->roworiented = flg; 1144 break; 1145 case MAT_KEEP_NONZERO_PATTERN: 1146 a->keepnonzeropattern = flg; 1147 break; 1148 case MAT_NEW_NONZERO_LOCATIONS: 1149 a->nonew = (flg ? 0 : 1); 1150 break; 1151 case MAT_NEW_NONZERO_LOCATION_ERR: 1152 a->nonew = (flg ? -1 : 0); 1153 break; 1154 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1155 a->nonew = (flg ? -2 : 0); 1156 break; 1157 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1158 a->nounused = (flg ? -1 : 0); 1159 break; 1160 case MAT_IGNORE_ZERO_ENTRIES: 1161 a->ignorezeroentries = flg; 1162 break; 1163 case MAT_SPD: 1164 case MAT_SYMMETRIC: 1165 case MAT_STRUCTURALLY_SYMMETRIC: 1166 case MAT_HERMITIAN: 1167 case MAT_SYMMETRY_ETERNAL: 1168 case MAT_STRUCTURE_ONLY: 1169 /* These options are handled directly by MatSetOption() */ 1170 break; 1171 case MAT_NEW_DIAGONALS: 1172 case MAT_IGNORE_OFF_PROC_ENTRIES: 1173 case MAT_USE_HASH_TABLE: 1174 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1175 break; 1176 case MAT_USE_INODES: 1177 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1178 break; 1179 case MAT_SUBMAT_SINGLEIS: 1180 A->submat_singleis = flg; 1181 break; 1182 default: 1183 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1184 } 1185 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1186 PetscFunctionReturn(0); 1187 } 1188 1189 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1190 { 1191 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1192 PetscErrorCode ierr; 1193 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1194 PetscScalar *aa=a->a,*x,zero=0.0; 1195 1196 PetscFunctionBegin; 1197 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1198 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1199 1200 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1201 PetscInt *diag=a->diag; 1202 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1203 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1204 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1205 PetscFunctionReturn(0); 1206 } 1207 1208 ierr = VecSet(v,zero);CHKERRQ(ierr); 1209 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1210 for (i=0; i<n; i++) { 1211 nz = ai[i+1] - ai[i]; 1212 if (!nz) x[i] = 0.0; 1213 for (j=ai[i]; j<ai[i+1]; j++) { 1214 if (aj[j] == i) { 1215 x[i] = aa[j]; 1216 break; 1217 } 1218 } 1219 } 1220 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1221 PetscFunctionReturn(0); 1222 } 1223 1224 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1225 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1226 { 1227 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1228 PetscScalar *y; 1229 const PetscScalar *x; 1230 PetscErrorCode ierr; 1231 PetscInt m = A->rmap->n; 1232 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1233 const MatScalar *v; 1234 PetscScalar alpha; 1235 PetscInt n,i,j; 1236 const PetscInt *idx,*ii,*ridx=NULL; 1237 Mat_CompressedRow cprow = a->compressedrow; 1238 PetscBool usecprow = cprow.use; 1239 #endif 1240 1241 PetscFunctionBegin; 1242 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1243 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1244 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1245 1246 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1247 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1248 #else 1249 if (usecprow) { 1250 m = cprow.nrows; 1251 ii = cprow.i; 1252 ridx = cprow.rindex; 1253 } else { 1254 ii = a->i; 1255 } 1256 for (i=0; i<m; i++) { 1257 idx = a->j + ii[i]; 1258 v = a->a + ii[i]; 1259 n = ii[i+1] - ii[i]; 1260 if (usecprow) { 1261 alpha = x[ridx[i]]; 1262 } else { 1263 alpha = x[i]; 1264 } 1265 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1266 } 1267 #endif 1268 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1269 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1270 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1271 PetscFunctionReturn(0); 1272 } 1273 1274 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1275 { 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1280 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1285 1286 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1287 { 1288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1289 PetscScalar *y; 1290 const PetscScalar *x; 1291 const MatScalar *aa; 1292 PetscErrorCode ierr; 1293 PetscInt m=A->rmap->n; 1294 const PetscInt *aj,*ii,*ridx=NULL; 1295 PetscInt n,i; 1296 PetscScalar sum; 1297 PetscBool usecprow=a->compressedrow.use; 1298 1299 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1300 #pragma disjoint(*x,*y,*aa) 1301 #endif 1302 1303 PetscFunctionBegin; 1304 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1305 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1306 ii = a->i; 1307 if (usecprow) { /* use compressed row format */ 1308 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1309 m = a->compressedrow.nrows; 1310 ii = a->compressedrow.i; 1311 ridx = a->compressedrow.rindex; 1312 for (i=0; i<m; i++) { 1313 n = ii[i+1] - ii[i]; 1314 aj = a->j + ii[i]; 1315 aa = a->a + ii[i]; 1316 sum = 0.0; 1317 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1318 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1319 y[*ridx++] = sum; 1320 } 1321 } else { /* do not use compressed row format */ 1322 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1323 aj = a->j; 1324 aa = a->a; 1325 fortranmultaij_(&m,x,ii,aj,aa,y); 1326 #else 1327 for (i=0; i<m; i++) { 1328 n = ii[i+1] - ii[i]; 1329 aj = a->j + ii[i]; 1330 aa = a->a + ii[i]; 1331 sum = 0.0; 1332 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1333 y[i] = sum; 1334 } 1335 #endif 1336 } 1337 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1338 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1339 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1344 { 1345 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1346 PetscScalar *y; 1347 const PetscScalar *x; 1348 const MatScalar *aa; 1349 PetscErrorCode ierr; 1350 PetscInt m=A->rmap->n; 1351 const PetscInt *aj,*ii,*ridx=NULL; 1352 PetscInt n,i,nonzerorow=0; 1353 PetscScalar sum; 1354 PetscBool usecprow=a->compressedrow.use; 1355 1356 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1357 #pragma disjoint(*x,*y,*aa) 1358 #endif 1359 1360 PetscFunctionBegin; 1361 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1362 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1363 if (usecprow) { /* use compressed row format */ 1364 m = a->compressedrow.nrows; 1365 ii = a->compressedrow.i; 1366 ridx = a->compressedrow.rindex; 1367 for (i=0; i<m; i++) { 1368 n = ii[i+1] - ii[i]; 1369 aj = a->j + ii[i]; 1370 aa = a->a + ii[i]; 1371 sum = 0.0; 1372 nonzerorow += (n>0); 1373 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1374 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1375 y[*ridx++] = sum; 1376 } 1377 } else { /* do not use compressed row format */ 1378 ii = a->i; 1379 for (i=0; i<m; i++) { 1380 n = ii[i+1] - ii[i]; 1381 aj = a->j + ii[i]; 1382 aa = a->a + ii[i]; 1383 sum = 0.0; 1384 nonzerorow += (n>0); 1385 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1386 y[i] = sum; 1387 } 1388 } 1389 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1390 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1391 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 1395 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1396 { 1397 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1398 PetscScalar *y,*z; 1399 const PetscScalar *x; 1400 const MatScalar *aa; 1401 PetscErrorCode ierr; 1402 PetscInt m = A->rmap->n,*aj,*ii; 1403 PetscInt n,i,*ridx=NULL; 1404 PetscScalar sum; 1405 PetscBool usecprow=a->compressedrow.use; 1406 1407 PetscFunctionBegin; 1408 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1409 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1410 if (usecprow) { /* use compressed row format */ 1411 if (zz != yy) { 1412 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1413 } 1414 m = a->compressedrow.nrows; 1415 ii = a->compressedrow.i; 1416 ridx = a->compressedrow.rindex; 1417 for (i=0; i<m; i++) { 1418 n = ii[i+1] - ii[i]; 1419 aj = a->j + ii[i]; 1420 aa = a->a + ii[i]; 1421 sum = y[*ridx]; 1422 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1423 z[*ridx++] = sum; 1424 } 1425 } else { /* do not use compressed row format */ 1426 ii = a->i; 1427 for (i=0; i<m; i++) { 1428 n = ii[i+1] - ii[i]; 1429 aj = a->j + ii[i]; 1430 aa = a->a + ii[i]; 1431 sum = y[i]; 1432 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1433 z[i] = sum; 1434 } 1435 } 1436 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1437 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1438 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1439 PetscFunctionReturn(0); 1440 } 1441 1442 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1443 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1444 { 1445 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1446 PetscScalar *y,*z; 1447 const PetscScalar *x; 1448 const MatScalar *aa; 1449 PetscErrorCode ierr; 1450 const PetscInt *aj,*ii,*ridx=NULL; 1451 PetscInt m = A->rmap->n,n,i; 1452 PetscScalar sum; 1453 PetscBool usecprow=a->compressedrow.use; 1454 1455 PetscFunctionBegin; 1456 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1457 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1458 if (usecprow) { /* use compressed row format */ 1459 if (zz != yy) { 1460 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1461 } 1462 m = a->compressedrow.nrows; 1463 ii = a->compressedrow.i; 1464 ridx = a->compressedrow.rindex; 1465 for (i=0; i<m; i++) { 1466 n = ii[i+1] - ii[i]; 1467 aj = a->j + ii[i]; 1468 aa = a->a + ii[i]; 1469 sum = y[*ridx]; 1470 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1471 z[*ridx++] = sum; 1472 } 1473 } else { /* do not use compressed row format */ 1474 ii = a->i; 1475 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1476 aj = a->j; 1477 aa = a->a; 1478 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1479 #else 1480 for (i=0; i<m; i++) { 1481 n = ii[i+1] - ii[i]; 1482 aj = a->j + ii[i]; 1483 aa = a->a + ii[i]; 1484 sum = y[i]; 1485 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1486 z[i] = sum; 1487 } 1488 #endif 1489 } 1490 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1491 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1492 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1493 PetscFunctionReturn(0); 1494 } 1495 1496 /* 1497 Adds diagonal pointers to sparse matrix structure. 1498 */ 1499 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1500 { 1501 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1502 PetscErrorCode ierr; 1503 PetscInt i,j,m = A->rmap->n; 1504 1505 PetscFunctionBegin; 1506 if (!a->diag) { 1507 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1508 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1509 } 1510 for (i=0; i<A->rmap->n; i++) { 1511 a->diag[i] = a->i[i+1]; 1512 for (j=a->i[i]; j<a->i[i+1]; j++) { 1513 if (a->j[j] == i) { 1514 a->diag[i] = j; 1515 break; 1516 } 1517 } 1518 } 1519 PetscFunctionReturn(0); 1520 } 1521 1522 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1523 { 1524 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1525 const PetscInt *diag = (const PetscInt*)a->diag; 1526 const PetscInt *ii = (const PetscInt*) a->i; 1527 PetscInt i,*mdiag = NULL; 1528 PetscErrorCode ierr; 1529 PetscInt cnt = 0; /* how many diagonals are missing */ 1530 1531 PetscFunctionBegin; 1532 if (!A->preallocated || !a->nz) { 1533 ierr = MatSeqAIJSetPreallocation(A,1,NULL);CHKERRQ(ierr); 1534 ierr = MatShift_Basic(A,v);CHKERRQ(ierr); 1535 PetscFunctionReturn(0); 1536 } 1537 1538 if (a->diagonaldense) { 1539 cnt = 0; 1540 } else { 1541 ierr = PetscCalloc1(A->rmap->n,&mdiag);CHKERRQ(ierr); 1542 for (i=0; i<A->rmap->n; i++) { 1543 if (diag[i] >= ii[i+1]) { 1544 cnt++; 1545 mdiag[i] = 1; 1546 } 1547 } 1548 } 1549 if (!cnt) { 1550 ierr = MatShift_Basic(A,v);CHKERRQ(ierr); 1551 } else { 1552 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1553 PetscInt *oldj = a->j, *oldi = a->i; 1554 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1555 1556 a->a = NULL; 1557 a->j = NULL; 1558 a->i = NULL; 1559 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1560 for (i=0; i<A->rmap->n; i++) { 1561 a->imax[i] += mdiag[i]; 1562 } 1563 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax);CHKERRQ(ierr); 1564 1565 /* copy old values into new matrix data structure */ 1566 for (i=0; i<A->rmap->n; i++) { 1567 ierr = MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES);CHKERRQ(ierr); 1568 ierr = MatSetValue(A,i,i,v,ADD_VALUES);CHKERRQ(ierr); 1569 } 1570 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1571 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1572 if (singlemalloc) { 1573 ierr = PetscFree3(olda,oldj,oldi);CHKERRQ(ierr); 1574 } else { 1575 if (free_a) {ierr = PetscFree(olda);CHKERRQ(ierr);} 1576 if (free_ij) {ierr = PetscFree(oldj);CHKERRQ(ierr);} 1577 if (free_ij) {ierr = PetscFree(oldi);CHKERRQ(ierr);} 1578 } 1579 } 1580 ierr = PetscFree(mdiag);CHKERRQ(ierr); 1581 a->diagonaldense = PETSC_TRUE; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 /* 1586 Checks for missing diagonals 1587 */ 1588 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1589 { 1590 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1591 PetscInt *diag,*ii = a->i,i; 1592 PetscErrorCode ierr; 1593 1594 PetscFunctionBegin; 1595 *missing = PETSC_FALSE; 1596 if (A->rmap->n > 0 && !ii) { 1597 *missing = PETSC_TRUE; 1598 if (d) *d = 0; 1599 ierr = PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");CHKERRQ(ierr); 1600 } else { 1601 diag = a->diag; 1602 for (i=0; i<A->rmap->n; i++) { 1603 if (diag[i] >= ii[i+1]) { 1604 *missing = PETSC_TRUE; 1605 if (d) *d = i; 1606 ierr = PetscInfo1(A,"Matrix is missing diagonal number %D\n",i);CHKERRQ(ierr); 1607 break; 1608 } 1609 } 1610 } 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #include <petscblaslapack.h> 1615 #include <petsc/private/kernels/blockinvert.h> 1616 1617 /* 1618 Note that values is allocated externally by the PC and then passed into this routine 1619 */ 1620 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1621 { 1622 PetscErrorCode ierr; 1623 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1624 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1625 const PetscReal shift = 0.0; 1626 PetscInt ipvt[5]; 1627 PetscScalar work[25],*v_work; 1628 1629 PetscFunctionBegin; 1630 allowzeropivot = PetscNot(A->erroriffailure); 1631 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1632 if (ncnt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %D doesn't match number matrix rows %D",ncnt,n); 1633 for (i=0; i<nblocks; i++) { 1634 bsizemax = PetscMax(bsizemax,bsizes[i]); 1635 } 1636 ierr = PetscMalloc1(bsizemax,&indx);CHKERRQ(ierr); 1637 if (bsizemax > 7) { 1638 ierr = PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots);CHKERRQ(ierr); 1639 } 1640 ncnt = 0; 1641 for (i=0; i<nblocks; i++) { 1642 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1643 ierr = MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag);CHKERRQ(ierr); 1644 switch (bsizes[i]) { 1645 case 1: 1646 *diag = 1.0/(*diag); 1647 break; 1648 case 2: 1649 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1650 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1651 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 1652 break; 1653 case 3: 1654 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1655 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1656 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 1657 break; 1658 case 4: 1659 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1660 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1661 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 1662 break; 1663 case 5: 1664 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1665 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1666 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 1667 break; 1668 case 6: 1669 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1670 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1671 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 1672 break; 1673 case 7: 1674 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1675 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1676 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 1677 break; 1678 default: 1679 ierr = PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1680 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1681 ierr = PetscKernel_A_gets_transpose_A_N(diag,bsizes[i]);CHKERRQ(ierr); 1682 } 1683 ncnt += bsizes[i]; 1684 diag += bsizes[i]*bsizes[i]; 1685 } 1686 if (bsizemax > 7) { 1687 ierr = PetscFree2(v_work,v_pivots);CHKERRQ(ierr); 1688 } 1689 ierr = PetscFree(indx);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 /* 1694 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1695 */ 1696 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1697 { 1698 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1699 PetscErrorCode ierr; 1700 PetscInt i,*diag,m = A->rmap->n; 1701 MatScalar *v = a->a; 1702 PetscScalar *idiag,*mdiag; 1703 1704 PetscFunctionBegin; 1705 if (a->idiagvalid) PetscFunctionReturn(0); 1706 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1707 diag = a->diag; 1708 if (!a->idiag) { 1709 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1710 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1711 v = a->a; 1712 } 1713 mdiag = a->mdiag; 1714 idiag = a->idiag; 1715 1716 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1717 for (i=0; i<m; i++) { 1718 mdiag[i] = v[diag[i]]; 1719 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1720 if (PetscRealPart(fshift)) { 1721 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1722 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1723 A->factorerror_zeropivot_value = 0.0; 1724 A->factorerror_zeropivot_row = i; 1725 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1726 } 1727 idiag[i] = 1.0/v[diag[i]]; 1728 } 1729 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1730 } else { 1731 for (i=0; i<m; i++) { 1732 mdiag[i] = v[diag[i]]; 1733 idiag[i] = omega/(fshift + v[diag[i]]); 1734 } 1735 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1736 } 1737 a->idiagvalid = PETSC_TRUE; 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1742 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1743 { 1744 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1745 PetscScalar *x,d,sum,*t,scale; 1746 const MatScalar *v,*idiag=0,*mdiag; 1747 const PetscScalar *b, *bs,*xb, *ts; 1748 PetscErrorCode ierr; 1749 PetscInt n,m = A->rmap->n,i; 1750 const PetscInt *idx,*diag; 1751 1752 PetscFunctionBegin; 1753 its = its*lits; 1754 1755 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1756 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1757 a->fshift = fshift; 1758 a->omega = omega; 1759 1760 diag = a->diag; 1761 t = a->ssor_work; 1762 idiag = a->idiag; 1763 mdiag = a->mdiag; 1764 1765 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1766 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1767 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1768 if (flag == SOR_APPLY_UPPER) { 1769 /* apply (U + D/omega) to the vector */ 1770 bs = b; 1771 for (i=0; i<m; i++) { 1772 d = fshift + mdiag[i]; 1773 n = a->i[i+1] - diag[i] - 1; 1774 idx = a->j + diag[i] + 1; 1775 v = a->a + diag[i] + 1; 1776 sum = b[i]*d/omega; 1777 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1778 x[i] = sum; 1779 } 1780 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1781 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1782 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1783 PetscFunctionReturn(0); 1784 } 1785 1786 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1787 else if (flag & SOR_EISENSTAT) { 1788 /* Let A = L + U + D; where L is lower trianglar, 1789 U is upper triangular, E = D/omega; This routine applies 1790 1791 (L + E)^{-1} A (U + E)^{-1} 1792 1793 to a vector efficiently using Eisenstat's trick. 1794 */ 1795 scale = (2.0/omega) - 1.0; 1796 1797 /* x = (E + U)^{-1} b */ 1798 for (i=m-1; i>=0; i--) { 1799 n = a->i[i+1] - diag[i] - 1; 1800 idx = a->j + diag[i] + 1; 1801 v = a->a + diag[i] + 1; 1802 sum = b[i]; 1803 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1804 x[i] = sum*idiag[i]; 1805 } 1806 1807 /* t = b - (2*E - D)x */ 1808 v = a->a; 1809 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1810 1811 /* t = (E + L)^{-1}t */ 1812 ts = t; 1813 diag = a->diag; 1814 for (i=0; i<m; i++) { 1815 n = diag[i] - a->i[i]; 1816 idx = a->j + a->i[i]; 1817 v = a->a + a->i[i]; 1818 sum = t[i]; 1819 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1820 t[i] = sum*idiag[i]; 1821 /* x = x + t */ 1822 x[i] += t[i]; 1823 } 1824 1825 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1826 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1827 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1828 PetscFunctionReturn(0); 1829 } 1830 if (flag & SOR_ZERO_INITIAL_GUESS) { 1831 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1832 for (i=0; i<m; i++) { 1833 n = diag[i] - a->i[i]; 1834 idx = a->j + a->i[i]; 1835 v = a->a + a->i[i]; 1836 sum = b[i]; 1837 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1838 t[i] = sum; 1839 x[i] = sum*idiag[i]; 1840 } 1841 xb = t; 1842 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1843 } else xb = b; 1844 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1845 for (i=m-1; i>=0; i--) { 1846 n = a->i[i+1] - diag[i] - 1; 1847 idx = a->j + diag[i] + 1; 1848 v = a->a + diag[i] + 1; 1849 sum = xb[i]; 1850 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1851 if (xb == b) { 1852 x[i] = sum*idiag[i]; 1853 } else { 1854 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1855 } 1856 } 1857 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1858 } 1859 its--; 1860 } 1861 while (its--) { 1862 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1863 for (i=0; i<m; i++) { 1864 /* lower */ 1865 n = diag[i] - a->i[i]; 1866 idx = a->j + a->i[i]; 1867 v = a->a + a->i[i]; 1868 sum = b[i]; 1869 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1870 t[i] = sum; /* save application of the lower-triangular part */ 1871 /* upper */ 1872 n = a->i[i+1] - diag[i] - 1; 1873 idx = a->j + diag[i] + 1; 1874 v = a->a + diag[i] + 1; 1875 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1876 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1877 } 1878 xb = t; 1879 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1880 } else xb = b; 1881 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1882 for (i=m-1; i>=0; i--) { 1883 sum = xb[i]; 1884 if (xb == b) { 1885 /* whole matrix (no checkpointing available) */ 1886 n = a->i[i+1] - a->i[i]; 1887 idx = a->j + a->i[i]; 1888 v = a->a + a->i[i]; 1889 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1890 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1891 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1892 n = a->i[i+1] - diag[i] - 1; 1893 idx = a->j + diag[i] + 1; 1894 v = a->a + diag[i] + 1; 1895 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1896 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1897 } 1898 } 1899 if (xb == b) { 1900 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1901 } else { 1902 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1903 } 1904 } 1905 } 1906 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1907 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 1912 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1913 { 1914 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1915 1916 PetscFunctionBegin; 1917 info->block_size = 1.0; 1918 info->nz_allocated = (double)a->maxnz; 1919 info->nz_used = (double)a->nz; 1920 info->nz_unneeded = (double)(a->maxnz - a->nz); 1921 info->assemblies = (double)A->num_ass; 1922 info->mallocs = (double)A->info.mallocs; 1923 info->memory = ((PetscObject)A)->mem; 1924 if (A->factortype) { 1925 info->fill_ratio_given = A->info.fill_ratio_given; 1926 info->fill_ratio_needed = A->info.fill_ratio_needed; 1927 info->factor_mallocs = A->info.factor_mallocs; 1928 } else { 1929 info->fill_ratio_given = 0; 1930 info->fill_ratio_needed = 0; 1931 info->factor_mallocs = 0; 1932 } 1933 PetscFunctionReturn(0); 1934 } 1935 1936 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1937 { 1938 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1939 PetscInt i,m = A->rmap->n - 1; 1940 PetscErrorCode ierr; 1941 const PetscScalar *xx; 1942 PetscScalar *bb; 1943 PetscInt d = 0; 1944 1945 PetscFunctionBegin; 1946 if (x && b) { 1947 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1948 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1949 for (i=0; i<N; i++) { 1950 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1951 bb[rows[i]] = diag*xx[rows[i]]; 1952 } 1953 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1954 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1955 } 1956 1957 if (a->keepnonzeropattern) { 1958 for (i=0; i<N; i++) { 1959 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1960 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1961 } 1962 if (diag != 0.0) { 1963 for (i=0; i<N; i++) { 1964 d = rows[i]; 1965 if (a->diag[d] >= a->i[d+1]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %D",d); 1966 } 1967 for (i=0; i<N; i++) { 1968 a->a[a->diag[rows[i]]] = diag; 1969 } 1970 } 1971 } else { 1972 if (diag != 0.0) { 1973 for (i=0; i<N; i++) { 1974 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1975 if (a->ilen[rows[i]] > 0) { 1976 a->ilen[rows[i]] = 1; 1977 a->a[a->i[rows[i]]] = diag; 1978 a->j[a->i[rows[i]]] = rows[i]; 1979 } else { /* in case row was completely empty */ 1980 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1981 } 1982 } 1983 } else { 1984 for (i=0; i<N; i++) { 1985 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1986 a->ilen[rows[i]] = 0; 1987 } 1988 } 1989 A->nonzerostate++; 1990 } 1991 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1992 PetscFunctionReturn(0); 1993 } 1994 1995 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1996 { 1997 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1998 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1999 PetscErrorCode ierr; 2000 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2001 const PetscScalar *xx; 2002 PetscScalar *bb; 2003 2004 PetscFunctionBegin; 2005 if (x && b) { 2006 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2007 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2008 vecs = PETSC_TRUE; 2009 } 2010 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2011 for (i=0; i<N; i++) { 2012 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2013 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2014 2015 zeroed[rows[i]] = PETSC_TRUE; 2016 } 2017 for (i=0; i<A->rmap->n; i++) { 2018 if (!zeroed[i]) { 2019 for (j=a->i[i]; j<a->i[i+1]; j++) { 2020 if (zeroed[a->j[j]]) { 2021 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2022 a->a[j] = 0.0; 2023 } 2024 } 2025 } else if (vecs) bb[i] = diag*xx[i]; 2026 } 2027 if (x && b) { 2028 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2029 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2030 } 2031 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2032 if (diag != 0.0) { 2033 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2034 if (missing) { 2035 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 2036 else { 2037 for (i=0; i<N; i++) { 2038 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 2039 } 2040 } 2041 } else { 2042 for (i=0; i<N; i++) { 2043 a->a[a->diag[rows[i]]] = diag; 2044 } 2045 } 2046 } 2047 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2052 { 2053 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2054 PetscInt *itmp; 2055 2056 PetscFunctionBegin; 2057 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2058 2059 *nz = a->i[row+1] - a->i[row]; 2060 if (v) *v = a->a + a->i[row]; 2061 if (idx) { 2062 itmp = a->j + a->i[row]; 2063 if (*nz) *idx = itmp; 2064 else *idx = 0; 2065 } 2066 PetscFunctionReturn(0); 2067 } 2068 2069 /* remove this function? */ 2070 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2071 { 2072 PetscFunctionBegin; 2073 PetscFunctionReturn(0); 2074 } 2075 2076 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2077 { 2078 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2079 MatScalar *v = a->a; 2080 PetscReal sum = 0.0; 2081 PetscErrorCode ierr; 2082 PetscInt i,j; 2083 2084 PetscFunctionBegin; 2085 if (type == NORM_FROBENIUS) { 2086 #if defined(PETSC_USE_REAL___FP16) 2087 PetscBLASInt one = 1,nz = a->nz; 2088 *nrm = BLASnrm2_(&nz,v,&one); 2089 #else 2090 for (i=0; i<a->nz; i++) { 2091 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2092 } 2093 *nrm = PetscSqrtReal(sum); 2094 #endif 2095 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2096 } else if (type == NORM_1) { 2097 PetscReal *tmp; 2098 PetscInt *jj = a->j; 2099 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2100 *nrm = 0.0; 2101 for (j=0; j<a->nz; j++) { 2102 tmp[*jj++] += PetscAbsScalar(*v); v++; 2103 } 2104 for (j=0; j<A->cmap->n; j++) { 2105 if (tmp[j] > *nrm) *nrm = tmp[j]; 2106 } 2107 ierr = PetscFree(tmp);CHKERRQ(ierr); 2108 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2109 } else if (type == NORM_INFINITY) { 2110 *nrm = 0.0; 2111 for (j=0; j<A->rmap->n; j++) { 2112 v = a->a + a->i[j]; 2113 sum = 0.0; 2114 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2115 sum += PetscAbsScalar(*v); v++; 2116 } 2117 if (sum > *nrm) *nrm = sum; 2118 } 2119 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2120 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2121 PetscFunctionReturn(0); 2122 } 2123 2124 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2125 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2126 { 2127 PetscErrorCode ierr; 2128 PetscInt i,j,anzj; 2129 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2130 PetscInt an=A->cmap->N,am=A->rmap->N; 2131 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2132 2133 PetscFunctionBegin; 2134 /* Allocate space for symbolic transpose info and work array */ 2135 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 2136 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2137 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2138 2139 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2140 /* Note: offset by 1 for fast conversion into csr format. */ 2141 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2142 /* Form ati for csr format of A^T. */ 2143 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2144 2145 /* Copy ati into atfill so we have locations of the next free space in atj */ 2146 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2147 2148 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2149 for (i=0;i<am;i++) { 2150 anzj = ai[i+1] - ai[i]; 2151 for (j=0;j<anzj;j++) { 2152 atj[atfill[*aj]] = i; 2153 atfill[*aj++] += 1; 2154 } 2155 } 2156 2157 /* Clean up temporary space and complete requests. */ 2158 ierr = PetscFree(atfill);CHKERRQ(ierr); 2159 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2160 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2161 2162 b = (Mat_SeqAIJ*)((*B)->data); 2163 b->free_a = PETSC_FALSE; 2164 b->free_ij = PETSC_TRUE; 2165 b->nonew = 0; 2166 PetscFunctionReturn(0); 2167 } 2168 2169 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2170 { 2171 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2172 Mat C; 2173 PetscErrorCode ierr; 2174 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2175 MatScalar *array = a->a; 2176 2177 PetscFunctionBegin; 2178 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 2179 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2180 2181 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2182 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2183 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2184 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2185 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2186 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2187 ierr = PetscFree(col);CHKERRQ(ierr); 2188 } else { 2189 C = *B; 2190 } 2191 2192 for (i=0; i<m; i++) { 2193 len = ai[i+1]-ai[i]; 2194 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2195 array += len; 2196 aj += len; 2197 } 2198 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2199 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2200 2201 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 2202 *B = C; 2203 } else { 2204 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2205 } 2206 PetscFunctionReturn(0); 2207 } 2208 2209 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2210 { 2211 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2212 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2213 MatScalar *va,*vb; 2214 PetscErrorCode ierr; 2215 PetscInt ma,na,mb,nb, i; 2216 2217 PetscFunctionBegin; 2218 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2219 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2220 if (ma!=nb || na!=mb) { 2221 *f = PETSC_FALSE; 2222 PetscFunctionReturn(0); 2223 } 2224 aii = aij->i; bii = bij->i; 2225 adx = aij->j; bdx = bij->j; 2226 va = aij->a; vb = bij->a; 2227 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2228 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2229 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2230 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2231 2232 *f = PETSC_TRUE; 2233 for (i=0; i<ma; i++) { 2234 while (aptr[i]<aii[i+1]) { 2235 PetscInt idc,idr; 2236 PetscScalar vc,vr; 2237 /* column/row index/value */ 2238 idc = adx[aptr[i]]; 2239 idr = bdx[bptr[idc]]; 2240 vc = va[aptr[i]]; 2241 vr = vb[bptr[idc]]; 2242 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2243 *f = PETSC_FALSE; 2244 goto done; 2245 } else { 2246 aptr[i]++; 2247 if (B || i!=idc) bptr[idc]++; 2248 } 2249 } 2250 } 2251 done: 2252 ierr = PetscFree(aptr);CHKERRQ(ierr); 2253 ierr = PetscFree(bptr);CHKERRQ(ierr); 2254 PetscFunctionReturn(0); 2255 } 2256 2257 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2258 { 2259 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2260 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2261 MatScalar *va,*vb; 2262 PetscErrorCode ierr; 2263 PetscInt ma,na,mb,nb, i; 2264 2265 PetscFunctionBegin; 2266 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2267 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2268 if (ma!=nb || na!=mb) { 2269 *f = PETSC_FALSE; 2270 PetscFunctionReturn(0); 2271 } 2272 aii = aij->i; bii = bij->i; 2273 adx = aij->j; bdx = bij->j; 2274 va = aij->a; vb = bij->a; 2275 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2276 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2277 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2278 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2279 2280 *f = PETSC_TRUE; 2281 for (i=0; i<ma; i++) { 2282 while (aptr[i]<aii[i+1]) { 2283 PetscInt idc,idr; 2284 PetscScalar vc,vr; 2285 /* column/row index/value */ 2286 idc = adx[aptr[i]]; 2287 idr = bdx[bptr[idc]]; 2288 vc = va[aptr[i]]; 2289 vr = vb[bptr[idc]]; 2290 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2291 *f = PETSC_FALSE; 2292 goto done; 2293 } else { 2294 aptr[i]++; 2295 if (B || i!=idc) bptr[idc]++; 2296 } 2297 } 2298 } 2299 done: 2300 ierr = PetscFree(aptr);CHKERRQ(ierr); 2301 ierr = PetscFree(bptr);CHKERRQ(ierr); 2302 PetscFunctionReturn(0); 2303 } 2304 2305 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2306 { 2307 PetscErrorCode ierr; 2308 2309 PetscFunctionBegin; 2310 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 2314 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2315 { 2316 PetscErrorCode ierr; 2317 2318 PetscFunctionBegin; 2319 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2324 { 2325 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2326 const PetscScalar *l,*r; 2327 PetscScalar x; 2328 MatScalar *v; 2329 PetscErrorCode ierr; 2330 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2331 const PetscInt *jj; 2332 2333 PetscFunctionBegin; 2334 if (ll) { 2335 /* The local size is used so that VecMPI can be passed to this routine 2336 by MatDiagonalScale_MPIAIJ */ 2337 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2338 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2339 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2340 v = a->a; 2341 for (i=0; i<m; i++) { 2342 x = l[i]; 2343 M = a->i[i+1] - a->i[i]; 2344 for (j=0; j<M; j++) (*v++) *= x; 2345 } 2346 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2347 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2348 } 2349 if (rr) { 2350 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2351 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2352 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2353 v = a->a; jj = a->j; 2354 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2355 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2356 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2357 } 2358 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2359 PetscFunctionReturn(0); 2360 } 2361 2362 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2363 { 2364 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2365 PetscErrorCode ierr; 2366 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2367 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2368 const PetscInt *irow,*icol; 2369 PetscInt nrows,ncols; 2370 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2371 MatScalar *a_new,*mat_a; 2372 Mat C; 2373 PetscBool stride; 2374 2375 PetscFunctionBegin; 2376 2377 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2378 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2379 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2380 2381 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2382 if (stride) { 2383 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2384 } else { 2385 first = 0; 2386 step = 0; 2387 } 2388 if (stride && step == 1) { 2389 /* special case of contiguous rows */ 2390 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2391 /* loop over new rows determining lens and starting points */ 2392 for (i=0; i<nrows; i++) { 2393 kstart = ai[irow[i]]; 2394 kend = kstart + ailen[irow[i]]; 2395 starts[i] = kstart; 2396 for (k=kstart; k<kend; k++) { 2397 if (aj[k] >= first) { 2398 starts[i] = k; 2399 break; 2400 } 2401 } 2402 sum = 0; 2403 while (k < kend) { 2404 if (aj[k++] >= first+ncols) break; 2405 sum++; 2406 } 2407 lens[i] = sum; 2408 } 2409 /* create submatrix */ 2410 if (scall == MAT_REUSE_MATRIX) { 2411 PetscInt n_cols,n_rows; 2412 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2413 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2414 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2415 C = *B; 2416 } else { 2417 PetscInt rbs,cbs; 2418 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2419 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2420 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2421 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2422 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2423 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2424 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2425 } 2426 c = (Mat_SeqAIJ*)C->data; 2427 2428 /* loop over rows inserting into submatrix */ 2429 a_new = c->a; 2430 j_new = c->j; 2431 i_new = c->i; 2432 2433 for (i=0; i<nrows; i++) { 2434 ii = starts[i]; 2435 lensi = lens[i]; 2436 for (k=0; k<lensi; k++) { 2437 *j_new++ = aj[ii+k] - first; 2438 } 2439 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2440 a_new += lensi; 2441 i_new[i+1] = i_new[i] + lensi; 2442 c->ilen[i] = lensi; 2443 } 2444 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2445 } else { 2446 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2447 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2448 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2449 for (i=0; i<ncols; i++) { 2450 #if defined(PETSC_USE_DEBUG) 2451 if (icol[i] >= oldcols) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%D] %D <= A->cmap->n %D",i,icol[i],oldcols); 2452 #endif 2453 smap[icol[i]] = i+1; 2454 } 2455 2456 /* determine lens of each row */ 2457 for (i=0; i<nrows; i++) { 2458 kstart = ai[irow[i]]; 2459 kend = kstart + a->ilen[irow[i]]; 2460 lens[i] = 0; 2461 for (k=kstart; k<kend; k++) { 2462 if (smap[aj[k]]) { 2463 lens[i]++; 2464 } 2465 } 2466 } 2467 /* Create and fill new matrix */ 2468 if (scall == MAT_REUSE_MATRIX) { 2469 PetscBool equal; 2470 2471 c = (Mat_SeqAIJ*)((*B)->data); 2472 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2473 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2474 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2475 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2476 C = *B; 2477 } else { 2478 PetscInt rbs,cbs; 2479 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2480 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2481 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2482 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2483 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2484 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2485 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2486 } 2487 c = (Mat_SeqAIJ*)(C->data); 2488 for (i=0; i<nrows; i++) { 2489 row = irow[i]; 2490 kstart = ai[row]; 2491 kend = kstart + a->ilen[row]; 2492 mat_i = c->i[i]; 2493 mat_j = c->j + mat_i; 2494 mat_a = c->a + mat_i; 2495 mat_ilen = c->ilen + i; 2496 for (k=kstart; k<kend; k++) { 2497 if ((tcol=smap[a->j[k]])) { 2498 *mat_j++ = tcol - 1; 2499 *mat_a++ = a->a[k]; 2500 (*mat_ilen)++; 2501 2502 } 2503 } 2504 } 2505 /* Free work space */ 2506 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2507 ierr = PetscFree(smap);CHKERRQ(ierr); 2508 ierr = PetscFree(lens);CHKERRQ(ierr); 2509 /* sort */ 2510 for (i = 0; i < nrows; i++) { 2511 PetscInt ilen; 2512 2513 mat_i = c->i[i]; 2514 mat_j = c->j + mat_i; 2515 mat_a = c->a + mat_i; 2516 ilen = c->ilen[i]; 2517 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2518 } 2519 } 2520 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2521 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2522 2523 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2524 *B = C; 2525 PetscFunctionReturn(0); 2526 } 2527 2528 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2529 { 2530 PetscErrorCode ierr; 2531 Mat B; 2532 2533 PetscFunctionBegin; 2534 if (scall == MAT_INITIAL_MATRIX) { 2535 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2536 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2537 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2538 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2539 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2540 *subMat = B; 2541 } else { 2542 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2543 } 2544 PetscFunctionReturn(0); 2545 } 2546 2547 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2548 { 2549 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2550 PetscErrorCode ierr; 2551 Mat outA; 2552 PetscBool row_identity,col_identity; 2553 2554 PetscFunctionBegin; 2555 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2556 2557 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2558 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2559 2560 outA = inA; 2561 outA->factortype = MAT_FACTOR_LU; 2562 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2563 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2564 2565 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2566 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2567 2568 a->row = row; 2569 2570 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2571 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2572 2573 a->col = col; 2574 2575 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2576 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2577 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2578 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2579 2580 if (!a->solve_work) { /* this matrix may have been factored before */ 2581 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2582 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2583 } 2584 2585 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2586 if (row_identity && col_identity) { 2587 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2588 } else { 2589 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2590 } 2591 PetscFunctionReturn(0); 2592 } 2593 2594 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2595 { 2596 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2597 PetscScalar oalpha = alpha; 2598 PetscErrorCode ierr; 2599 PetscBLASInt one = 1,bnz; 2600 2601 PetscFunctionBegin; 2602 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2603 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2604 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2605 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2606 PetscFunctionReturn(0); 2607 } 2608 2609 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2610 { 2611 PetscErrorCode ierr; 2612 PetscInt i; 2613 2614 PetscFunctionBegin; 2615 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2616 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2617 2618 for (i=0; i<submatj->nrqr; ++i) { 2619 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2620 } 2621 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2622 2623 if (submatj->rbuf1) { 2624 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2625 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2626 } 2627 2628 for (i=0; i<submatj->nrqs; ++i) { 2629 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2630 } 2631 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2632 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2633 } 2634 2635 #if defined(PETSC_USE_CTABLE) 2636 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2637 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2638 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2639 #else 2640 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2641 #endif 2642 2643 if (!submatj->allcolumns) { 2644 #if defined(PETSC_USE_CTABLE) 2645 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2646 #else 2647 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2648 #endif 2649 } 2650 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2651 2652 ierr = PetscFree(submatj);CHKERRQ(ierr); 2653 PetscFunctionReturn(0); 2654 } 2655 2656 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2657 { 2658 PetscErrorCode ierr; 2659 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2660 Mat_SubSppt *submatj = c->submatis1; 2661 2662 PetscFunctionBegin; 2663 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2664 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2665 PetscFunctionReturn(0); 2666 } 2667 2668 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2669 { 2670 PetscErrorCode ierr; 2671 PetscInt i; 2672 Mat C; 2673 Mat_SeqAIJ *c; 2674 Mat_SubSppt *submatj; 2675 2676 PetscFunctionBegin; 2677 for (i=0; i<n; i++) { 2678 C = (*mat)[i]; 2679 c = (Mat_SeqAIJ*)C->data; 2680 submatj = c->submatis1; 2681 if (submatj) { 2682 if (--((PetscObject)C)->refct <= 0) { 2683 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2684 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2685 ierr = PetscFree(C->defaultvectype);CHKERRQ(ierr); 2686 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2687 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2688 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2689 } 2690 } else { 2691 ierr = MatDestroy(&C);CHKERRQ(ierr); 2692 } 2693 } 2694 2695 /* Destroy Dummy submatrices created for reuse */ 2696 ierr = MatDestroySubMatrices_Dummy(n,mat);CHKERRQ(ierr); 2697 2698 ierr = PetscFree(*mat);CHKERRQ(ierr); 2699 PetscFunctionReturn(0); 2700 } 2701 2702 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2703 { 2704 PetscErrorCode ierr; 2705 PetscInt i; 2706 2707 PetscFunctionBegin; 2708 if (scall == MAT_INITIAL_MATRIX) { 2709 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2710 } 2711 2712 for (i=0; i<n; i++) { 2713 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2714 } 2715 PetscFunctionReturn(0); 2716 } 2717 2718 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2719 { 2720 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2721 PetscErrorCode ierr; 2722 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2723 const PetscInt *idx; 2724 PetscInt start,end,*ai,*aj; 2725 PetscBT table; 2726 2727 PetscFunctionBegin; 2728 m = A->rmap->n; 2729 ai = a->i; 2730 aj = a->j; 2731 2732 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2733 2734 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2735 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2736 2737 for (i=0; i<is_max; i++) { 2738 /* Initialize the two local arrays */ 2739 isz = 0; 2740 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2741 2742 /* Extract the indices, assume there can be duplicate entries */ 2743 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2744 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2745 2746 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2747 for (j=0; j<n; ++j) { 2748 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2749 } 2750 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2751 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2752 2753 k = 0; 2754 for (j=0; j<ov; j++) { /* for each overlap */ 2755 n = isz; 2756 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2757 row = nidx[k]; 2758 start = ai[row]; 2759 end = ai[row+1]; 2760 for (l = start; l<end; l++) { 2761 val = aj[l]; 2762 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2763 } 2764 } 2765 } 2766 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2767 } 2768 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2769 ierr = PetscFree(nidx);CHKERRQ(ierr); 2770 PetscFunctionReturn(0); 2771 } 2772 2773 /* -------------------------------------------------------------- */ 2774 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2775 { 2776 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2777 PetscErrorCode ierr; 2778 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2779 const PetscInt *row,*col; 2780 PetscInt *cnew,j,*lens; 2781 IS icolp,irowp; 2782 PetscInt *cwork = NULL; 2783 PetscScalar *vwork = NULL; 2784 2785 PetscFunctionBegin; 2786 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2787 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2788 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2789 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2790 2791 /* determine lengths of permuted rows */ 2792 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2793 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2794 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2795 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2796 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2797 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2798 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2799 ierr = PetscFree(lens);CHKERRQ(ierr); 2800 2801 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2802 for (i=0; i<m; i++) { 2803 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2804 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2805 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2806 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2807 } 2808 ierr = PetscFree(cnew);CHKERRQ(ierr); 2809 2810 (*B)->assembled = PETSC_FALSE; 2811 2812 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2813 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2814 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2815 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2816 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2817 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2818 PetscFunctionReturn(0); 2819 } 2820 2821 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2822 { 2823 PetscErrorCode ierr; 2824 2825 PetscFunctionBegin; 2826 /* If the two matrices have the same copy implementation, use fast copy. */ 2827 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2828 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2829 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2830 2831 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"); 2832 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2833 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2834 } else { 2835 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2836 } 2837 PetscFunctionReturn(0); 2838 } 2839 2840 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2841 { 2842 PetscErrorCode ierr; 2843 2844 PetscFunctionBegin; 2845 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2846 PetscFunctionReturn(0); 2847 } 2848 2849 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2850 { 2851 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2852 2853 PetscFunctionBegin; 2854 *array = a->a; 2855 PetscFunctionReturn(0); 2856 } 2857 2858 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2859 { 2860 PetscFunctionBegin; 2861 PetscFunctionReturn(0); 2862 } 2863 2864 /* 2865 Computes the number of nonzeros per row needed for preallocation when X and Y 2866 have different nonzero structure. 2867 */ 2868 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2869 { 2870 PetscInt i,j,k,nzx,nzy; 2871 2872 PetscFunctionBegin; 2873 /* Set the number of nonzeros in the new matrix */ 2874 for (i=0; i<m; i++) { 2875 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2876 nzx = xi[i+1] - xi[i]; 2877 nzy = yi[i+1] - yi[i]; 2878 nnz[i] = 0; 2879 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2880 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2881 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2882 nnz[i]++; 2883 } 2884 for (; k<nzy; k++) nnz[i]++; 2885 } 2886 PetscFunctionReturn(0); 2887 } 2888 2889 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2890 { 2891 PetscInt m = Y->rmap->N; 2892 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2893 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2894 PetscErrorCode ierr; 2895 2896 PetscFunctionBegin; 2897 /* Set the number of nonzeros in the new matrix */ 2898 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2899 PetscFunctionReturn(0); 2900 } 2901 2902 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2903 { 2904 PetscErrorCode ierr; 2905 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2906 PetscBLASInt one=1,bnz; 2907 2908 PetscFunctionBegin; 2909 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2910 if (str == SAME_NONZERO_PATTERN) { 2911 PetscScalar alpha = a; 2912 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2913 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2914 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2915 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2916 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2917 } else { 2918 Mat B; 2919 PetscInt *nnz; 2920 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2921 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2922 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2923 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2924 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2925 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2926 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2927 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2928 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2929 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2930 ierr = PetscFree(nnz);CHKERRQ(ierr); 2931 } 2932 PetscFunctionReturn(0); 2933 } 2934 2935 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2936 { 2937 #if defined(PETSC_USE_COMPLEX) 2938 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2939 PetscInt i,nz; 2940 PetscScalar *a; 2941 2942 PetscFunctionBegin; 2943 nz = aij->nz; 2944 a = aij->a; 2945 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2946 #else 2947 PetscFunctionBegin; 2948 #endif 2949 PetscFunctionReturn(0); 2950 } 2951 2952 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2953 { 2954 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2955 PetscErrorCode ierr; 2956 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2957 PetscReal atmp; 2958 PetscScalar *x; 2959 MatScalar *aa; 2960 2961 PetscFunctionBegin; 2962 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2963 aa = a->a; 2964 ai = a->i; 2965 aj = a->j; 2966 2967 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2968 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2969 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2970 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2971 for (i=0; i<m; i++) { 2972 ncols = ai[1] - ai[0]; ai++; 2973 x[i] = 0.0; 2974 for (j=0; j<ncols; j++) { 2975 atmp = PetscAbsScalar(*aa); 2976 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2977 aa++; aj++; 2978 } 2979 } 2980 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2981 PetscFunctionReturn(0); 2982 } 2983 2984 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2985 { 2986 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2987 PetscErrorCode ierr; 2988 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2989 PetscScalar *x; 2990 MatScalar *aa; 2991 2992 PetscFunctionBegin; 2993 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2994 aa = a->a; 2995 ai = a->i; 2996 aj = a->j; 2997 2998 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2999 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3000 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3001 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3002 for (i=0; i<m; i++) { 3003 ncols = ai[1] - ai[0]; ai++; 3004 if (ncols == A->cmap->n) { /* row is dense */ 3005 x[i] = *aa; if (idx) idx[i] = 0; 3006 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3007 x[i] = 0.0; 3008 if (idx) { 3009 idx[i] = 0; /* in case ncols is zero */ 3010 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 3011 if (aj[j] > j) { 3012 idx[i] = j; 3013 break; 3014 } 3015 } 3016 } 3017 } 3018 for (j=0; j<ncols; j++) { 3019 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3020 aa++; aj++; 3021 } 3022 } 3023 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3024 PetscFunctionReturn(0); 3025 } 3026 3027 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3028 { 3029 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3030 PetscErrorCode ierr; 3031 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3032 PetscReal atmp; 3033 PetscScalar *x; 3034 MatScalar *aa; 3035 3036 PetscFunctionBegin; 3037 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3038 aa = a->a; 3039 ai = a->i; 3040 aj = a->j; 3041 3042 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3043 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3044 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3045 if (n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %D vs. %D rows", A->rmap->n, n); 3046 for (i=0; i<m; i++) { 3047 ncols = ai[1] - ai[0]; ai++; 3048 if (ncols) { 3049 /* Get first nonzero */ 3050 for (j = 0; j < ncols; j++) { 3051 atmp = PetscAbsScalar(aa[j]); 3052 if (atmp > 1.0e-12) { 3053 x[i] = atmp; 3054 if (idx) idx[i] = aj[j]; 3055 break; 3056 } 3057 } 3058 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3059 } else { 3060 x[i] = 0.0; if (idx) idx[i] = 0; 3061 } 3062 for (j = 0; j < ncols; j++) { 3063 atmp = PetscAbsScalar(*aa); 3064 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3065 aa++; aj++; 3066 } 3067 } 3068 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3069 PetscFunctionReturn(0); 3070 } 3071 3072 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3073 { 3074 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3075 PetscErrorCode ierr; 3076 PetscInt i,j,m = A->rmap->n,ncols,n; 3077 const PetscInt *ai,*aj; 3078 PetscScalar *x; 3079 const MatScalar *aa; 3080 3081 PetscFunctionBegin; 3082 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3083 aa = a->a; 3084 ai = a->i; 3085 aj = a->j; 3086 3087 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3088 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3089 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3090 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3091 for (i=0; i<m; i++) { 3092 ncols = ai[1] - ai[0]; ai++; 3093 if (ncols == A->cmap->n) { /* row is dense */ 3094 x[i] = *aa; if (idx) idx[i] = 0; 3095 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3096 x[i] = 0.0; 3097 if (idx) { /* find first implicit 0.0 in the row */ 3098 idx[i] = 0; /* in case ncols is zero */ 3099 for (j=0; j<ncols; j++) { 3100 if (aj[j] > j) { 3101 idx[i] = j; 3102 break; 3103 } 3104 } 3105 } 3106 } 3107 for (j=0; j<ncols; j++) { 3108 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3109 aa++; aj++; 3110 } 3111 } 3112 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3113 PetscFunctionReturn(0); 3114 } 3115 3116 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3117 { 3118 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3119 PetscErrorCode ierr; 3120 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3121 MatScalar *diag,work[25],*v_work; 3122 const PetscReal shift = 0.0; 3123 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3124 3125 PetscFunctionBegin; 3126 allowzeropivot = PetscNot(A->erroriffailure); 3127 if (a->ibdiagvalid) { 3128 if (values) *values = a->ibdiag; 3129 PetscFunctionReturn(0); 3130 } 3131 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3132 if (!a->ibdiag) { 3133 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3134 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3135 } 3136 diag = a->ibdiag; 3137 if (values) *values = a->ibdiag; 3138 /* factor and invert each block */ 3139 switch (bs) { 3140 case 1: 3141 for (i=0; i<mbs; i++) { 3142 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3143 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3144 if (allowzeropivot) { 3145 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3146 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3147 A->factorerror_zeropivot_row = i; 3148 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3149 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3150 } 3151 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3152 } 3153 break; 3154 case 2: 3155 for (i=0; i<mbs; i++) { 3156 ij[0] = 2*i; ij[1] = 2*i + 1; 3157 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3158 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3159 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3160 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3161 diag += 4; 3162 } 3163 break; 3164 case 3: 3165 for (i=0; i<mbs; i++) { 3166 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3167 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3168 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3169 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3170 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3171 diag += 9; 3172 } 3173 break; 3174 case 4: 3175 for (i=0; i<mbs; i++) { 3176 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3177 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3178 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3179 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3180 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3181 diag += 16; 3182 } 3183 break; 3184 case 5: 3185 for (i=0; i<mbs; i++) { 3186 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3187 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3188 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3189 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3190 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3191 diag += 25; 3192 } 3193 break; 3194 case 6: 3195 for (i=0; i<mbs; i++) { 3196 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3197 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3198 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3199 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3200 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3201 diag += 36; 3202 } 3203 break; 3204 case 7: 3205 for (i=0; i<mbs; i++) { 3206 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3207 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3208 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3209 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3210 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3211 diag += 49; 3212 } 3213 break; 3214 default: 3215 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3216 for (i=0; i<mbs; i++) { 3217 for (j=0; j<bs; j++) { 3218 IJ[j] = bs*i + j; 3219 } 3220 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3221 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3222 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3223 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3224 diag += bs2; 3225 } 3226 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3227 } 3228 a->ibdiagvalid = PETSC_TRUE; 3229 PetscFunctionReturn(0); 3230 } 3231 3232 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3233 { 3234 PetscErrorCode ierr; 3235 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3236 PetscScalar a; 3237 PetscInt m,n,i,j,col; 3238 3239 PetscFunctionBegin; 3240 if (!x->assembled) { 3241 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3242 for (i=0; i<m; i++) { 3243 for (j=0; j<aij->imax[i]; j++) { 3244 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3245 col = (PetscInt)(n*PetscRealPart(a)); 3246 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3247 } 3248 } 3249 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3250 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3251 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3252 PetscFunctionReturn(0); 3253 } 3254 3255 /* -------------------------------------------------------------------*/ 3256 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3257 MatGetRow_SeqAIJ, 3258 MatRestoreRow_SeqAIJ, 3259 MatMult_SeqAIJ, 3260 /* 4*/ MatMultAdd_SeqAIJ, 3261 MatMultTranspose_SeqAIJ, 3262 MatMultTransposeAdd_SeqAIJ, 3263 0, 3264 0, 3265 0, 3266 /* 10*/ 0, 3267 MatLUFactor_SeqAIJ, 3268 0, 3269 MatSOR_SeqAIJ, 3270 MatTranspose_SeqAIJ_FAST, 3271 /*1 5*/ MatGetInfo_SeqAIJ, 3272 MatEqual_SeqAIJ, 3273 MatGetDiagonal_SeqAIJ, 3274 MatDiagonalScale_SeqAIJ, 3275 MatNorm_SeqAIJ, 3276 /* 20*/ 0, 3277 MatAssemblyEnd_SeqAIJ, 3278 MatSetOption_SeqAIJ, 3279 MatZeroEntries_SeqAIJ, 3280 /* 24*/ MatZeroRows_SeqAIJ, 3281 0, 3282 0, 3283 0, 3284 0, 3285 /* 29*/ MatSetUp_SeqAIJ, 3286 0, 3287 0, 3288 0, 3289 0, 3290 /* 34*/ MatDuplicate_SeqAIJ, 3291 0, 3292 0, 3293 MatILUFactor_SeqAIJ, 3294 0, 3295 /* 39*/ MatAXPY_SeqAIJ, 3296 MatCreateSubMatrices_SeqAIJ, 3297 MatIncreaseOverlap_SeqAIJ, 3298 MatGetValues_SeqAIJ, 3299 MatCopy_SeqAIJ, 3300 /* 44*/ MatGetRowMax_SeqAIJ, 3301 MatScale_SeqAIJ, 3302 MatShift_SeqAIJ, 3303 MatDiagonalSet_SeqAIJ, 3304 MatZeroRowsColumns_SeqAIJ, 3305 /* 49*/ MatSetRandom_SeqAIJ, 3306 MatGetRowIJ_SeqAIJ, 3307 MatRestoreRowIJ_SeqAIJ, 3308 MatGetColumnIJ_SeqAIJ, 3309 MatRestoreColumnIJ_SeqAIJ, 3310 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3311 0, 3312 0, 3313 MatPermute_SeqAIJ, 3314 0, 3315 /* 59*/ 0, 3316 MatDestroy_SeqAIJ, 3317 MatView_SeqAIJ, 3318 0, 3319 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3320 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3321 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3322 0, 3323 0, 3324 0, 3325 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3326 MatGetRowMinAbs_SeqAIJ, 3327 0, 3328 0, 3329 0, 3330 /* 74*/ 0, 3331 MatFDColoringApply_AIJ, 3332 0, 3333 0, 3334 0, 3335 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3336 0, 3337 0, 3338 0, 3339 MatLoad_SeqAIJ, 3340 /* 84*/ MatIsSymmetric_SeqAIJ, 3341 MatIsHermitian_SeqAIJ, 3342 0, 3343 0, 3344 0, 3345 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3346 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3347 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3348 MatPtAP_SeqAIJ_SeqAIJ, 3349 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3350 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3351 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3352 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3353 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3354 0, 3355 /* 99*/ 0, 3356 0, 3357 0, 3358 MatConjugate_SeqAIJ, 3359 0, 3360 /*104*/ MatSetValuesRow_SeqAIJ, 3361 MatRealPart_SeqAIJ, 3362 MatImaginaryPart_SeqAIJ, 3363 0, 3364 0, 3365 /*109*/ MatMatSolve_SeqAIJ, 3366 0, 3367 MatGetRowMin_SeqAIJ, 3368 0, 3369 MatMissingDiagonal_SeqAIJ, 3370 /*114*/ 0, 3371 0, 3372 0, 3373 0, 3374 0, 3375 /*119*/ 0, 3376 0, 3377 0, 3378 0, 3379 MatGetMultiProcBlock_SeqAIJ, 3380 /*124*/ MatFindNonzeroRows_SeqAIJ, 3381 MatGetColumnNorms_SeqAIJ, 3382 MatInvertBlockDiagonal_SeqAIJ, 3383 MatInvertVariableBlockDiagonal_SeqAIJ, 3384 0, 3385 /*129*/ 0, 3386 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3387 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3388 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3389 MatTransposeColoringCreate_SeqAIJ, 3390 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3391 MatTransColoringApplyDenToSp_SeqAIJ, 3392 MatRARt_SeqAIJ_SeqAIJ, 3393 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3394 MatRARtNumeric_SeqAIJ_SeqAIJ, 3395 /*139*/0, 3396 0, 3397 0, 3398 MatFDColoringSetUp_SeqXAIJ, 3399 MatFindOffBlockDiagonalEntries_SeqAIJ, 3400 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3401 MatDestroySubMatrices_SeqAIJ 3402 }; 3403 3404 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3405 { 3406 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3407 PetscInt i,nz,n; 3408 3409 PetscFunctionBegin; 3410 nz = aij->maxnz; 3411 n = mat->rmap->n; 3412 for (i=0; i<nz; i++) { 3413 aij->j[i] = indices[i]; 3414 } 3415 aij->nz = nz; 3416 for (i=0; i<n; i++) { 3417 aij->ilen[i] = aij->imax[i]; 3418 } 3419 PetscFunctionReturn(0); 3420 } 3421 3422 /* 3423 * When a sparse matrix has many zero colums, we should compact them out to save the space 3424 * This happens in MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3425 * */ 3426 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3427 { 3428 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3429 #if defined(PETSC_USE_CTABLE) 3430 PetscTable gid1_lid1; 3431 PetscTablePosition tpos; 3432 PetscInt gid,lid,i,j,ncols,ec; 3433 PetscInt *garray; 3434 #endif 3435 PetscErrorCode ierr; 3436 3437 PetscFunctionBegin; 3438 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3439 PetscValidPointer(mapping,2); 3440 #if defined(PETSC_USE_CTABLE) 3441 /* use a table */ 3442 ierr = PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr); 3443 ec = 0; 3444 for (i=0; i<mat->rmap->n; i++) { 3445 ncols = aij->i[i+1] - aij->i[i]; 3446 for (j=0; j<ncols; j++) { 3447 PetscInt data,gid1 = aij->j[aij->i[i] + j] + 1; 3448 ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr); 3449 if (!data) { 3450 /* one based table */ 3451 ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr); 3452 } 3453 } 3454 } 3455 /* form array of columns we need */ 3456 ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr); 3457 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 3458 while (tpos) { 3459 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr); 3460 gid--; 3461 lid--; 3462 garray[lid] = gid; 3463 } 3464 ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */ 3465 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 3466 for (i=0; i<ec; i++) { 3467 ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 3468 } 3469 /* compact out the extra columns in B */ 3470 for (i=0; i<mat->rmap->n; i++) { 3471 ncols = aij->i[i+1] - aij->i[i]; 3472 for (j=0; j<ncols; j++) { 3473 PetscInt gid1 = aij->j[aij->i[i] + j] + 1; 3474 ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr); 3475 lid--; 3476 aij->j[aij->i[i] + j] = lid; 3477 } 3478 } 3479 mat->cmap->n = mat->cmap->N = ec; 3480 mat->cmap->bs = 1; 3481 3482 ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr); 3483 ierr = PetscLayoutSetUp((mat->cmap));CHKERRQ(ierr); 3484 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);CHKERRQ(ierr); 3485 ierr = ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);CHKERRQ(ierr); 3486 #else 3487 SETERRQ(PETSC_COMM_SELF,"Require ctable \n"); 3488 #endif 3489 PetscFunctionReturn(0); 3490 } 3491 3492 /*@ 3493 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3494 in the matrix. 3495 3496 Input Parameters: 3497 + mat - the SeqAIJ matrix 3498 - indices - the column indices 3499 3500 Level: advanced 3501 3502 Notes: 3503 This can be called if you have precomputed the nonzero structure of the 3504 matrix and want to provide it to the matrix object to improve the performance 3505 of the MatSetValues() operation. 3506 3507 You MUST have set the correct numbers of nonzeros per row in the call to 3508 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3509 3510 MUST be called before any calls to MatSetValues(); 3511 3512 The indices should start with zero, not one. 3513 3514 @*/ 3515 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3516 { 3517 PetscErrorCode ierr; 3518 3519 PetscFunctionBegin; 3520 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3521 PetscValidPointer(indices,2); 3522 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3523 PetscFunctionReturn(0); 3524 } 3525 3526 /* ----------------------------------------------------------------------------------------*/ 3527 3528 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3529 { 3530 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3531 PetscErrorCode ierr; 3532 size_t nz = aij->i[mat->rmap->n]; 3533 3534 PetscFunctionBegin; 3535 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3536 3537 /* allocate space for values if not already there */ 3538 if (!aij->saved_values) { 3539 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3540 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3541 } 3542 3543 /* copy values over */ 3544 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3545 PetscFunctionReturn(0); 3546 } 3547 3548 /*@ 3549 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3550 example, reuse of the linear part of a Jacobian, while recomputing the 3551 nonlinear portion. 3552 3553 Collect on Mat 3554 3555 Input Parameters: 3556 . mat - the matrix (currently only AIJ matrices support this option) 3557 3558 Level: advanced 3559 3560 Common Usage, with SNESSolve(): 3561 $ Create Jacobian matrix 3562 $ Set linear terms into matrix 3563 $ Apply boundary conditions to matrix, at this time matrix must have 3564 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3565 $ boundary conditions again will not change the nonzero structure 3566 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3567 $ ierr = MatStoreValues(mat); 3568 $ Call SNESSetJacobian() with matrix 3569 $ In your Jacobian routine 3570 $ ierr = MatRetrieveValues(mat); 3571 $ Set nonlinear terms in matrix 3572 3573 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3574 $ // build linear portion of Jacobian 3575 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3576 $ ierr = MatStoreValues(mat); 3577 $ loop over nonlinear iterations 3578 $ ierr = MatRetrieveValues(mat); 3579 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3580 $ // call MatAssemblyBegin/End() on matrix 3581 $ Solve linear system with Jacobian 3582 $ endloop 3583 3584 Notes: 3585 Matrix must already be assemblied before calling this routine 3586 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3587 calling this routine. 3588 3589 When this is called multiple times it overwrites the previous set of stored values 3590 and does not allocated additional space. 3591 3592 .seealso: MatRetrieveValues() 3593 3594 @*/ 3595 PetscErrorCode MatStoreValues(Mat mat) 3596 { 3597 PetscErrorCode ierr; 3598 3599 PetscFunctionBegin; 3600 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3601 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3602 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3603 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3604 PetscFunctionReturn(0); 3605 } 3606 3607 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3608 { 3609 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3610 PetscErrorCode ierr; 3611 PetscInt nz = aij->i[mat->rmap->n]; 3612 3613 PetscFunctionBegin; 3614 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3615 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3616 /* copy values over */ 3617 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3618 PetscFunctionReturn(0); 3619 } 3620 3621 /*@ 3622 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3623 example, reuse of the linear part of a Jacobian, while recomputing the 3624 nonlinear portion. 3625 3626 Collect on Mat 3627 3628 Input Parameters: 3629 . mat - the matrix (currently only AIJ matrices support this option) 3630 3631 Level: advanced 3632 3633 .seealso: MatStoreValues() 3634 3635 @*/ 3636 PetscErrorCode MatRetrieveValues(Mat mat) 3637 { 3638 PetscErrorCode ierr; 3639 3640 PetscFunctionBegin; 3641 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3642 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3643 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3644 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3645 PetscFunctionReturn(0); 3646 } 3647 3648 3649 /* --------------------------------------------------------------------------------*/ 3650 /*@C 3651 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3652 (the default parallel PETSc format). For good matrix assembly performance 3653 the user should preallocate the matrix storage by setting the parameter nz 3654 (or the array nnz). By setting these parameters accurately, performance 3655 during matrix assembly can be increased by more than a factor of 50. 3656 3657 Collective on MPI_Comm 3658 3659 Input Parameters: 3660 + comm - MPI communicator, set to PETSC_COMM_SELF 3661 . m - number of rows 3662 . n - number of columns 3663 . nz - number of nonzeros per row (same for all rows) 3664 - nnz - array containing the number of nonzeros in the various rows 3665 (possibly different for each row) or NULL 3666 3667 Output Parameter: 3668 . A - the matrix 3669 3670 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3671 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3672 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3673 3674 Notes: 3675 If nnz is given then nz is ignored 3676 3677 The AIJ format (also called the Yale sparse matrix format or 3678 compressed row storage), is fully compatible with standard Fortran 77 3679 storage. That is, the stored row and column indices can begin at 3680 either one (as in Fortran) or zero. See the users' manual for details. 3681 3682 Specify the preallocated storage with either nz or nnz (not both). 3683 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3684 allocation. For large problems you MUST preallocate memory or you 3685 will get TERRIBLE performance, see the users' manual chapter on matrices. 3686 3687 By default, this format uses inodes (identical nodes) when possible, to 3688 improve numerical efficiency of matrix-vector products and solves. We 3689 search for consecutive rows with the same nonzero structure, thereby 3690 reusing matrix information to achieve increased efficiency. 3691 3692 Options Database Keys: 3693 + -mat_no_inode - Do not use inodes 3694 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3695 3696 Level: intermediate 3697 3698 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3699 3700 @*/ 3701 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3702 { 3703 PetscErrorCode ierr; 3704 3705 PetscFunctionBegin; 3706 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3707 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3708 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3709 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3710 PetscFunctionReturn(0); 3711 } 3712 3713 /*@C 3714 MatSeqAIJSetPreallocation - For good matrix assembly performance 3715 the user should preallocate the matrix storage by setting the parameter nz 3716 (or the array nnz). By setting these parameters accurately, performance 3717 during matrix assembly can be increased by more than a factor of 50. 3718 3719 Collective on MPI_Comm 3720 3721 Input Parameters: 3722 + B - The matrix 3723 . nz - number of nonzeros per row (same for all rows) 3724 - nnz - array containing the number of nonzeros in the various rows 3725 (possibly different for each row) or NULL 3726 3727 Notes: 3728 If nnz is given then nz is ignored 3729 3730 The AIJ format (also called the Yale sparse matrix format or 3731 compressed row storage), is fully compatible with standard Fortran 77 3732 storage. That is, the stored row and column indices can begin at 3733 either one (as in Fortran) or zero. See the users' manual for details. 3734 3735 Specify the preallocated storage with either nz or nnz (not both). 3736 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3737 allocation. For large problems you MUST preallocate memory or you 3738 will get TERRIBLE performance, see the users' manual chapter on matrices. 3739 3740 You can call MatGetInfo() to get information on how effective the preallocation was; 3741 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3742 You can also run with the option -info and look for messages with the string 3743 malloc in them to see if additional memory allocation was needed. 3744 3745 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3746 entries or columns indices 3747 3748 By default, this format uses inodes (identical nodes) when possible, to 3749 improve numerical efficiency of matrix-vector products and solves. We 3750 search for consecutive rows with the same nonzero structure, thereby 3751 reusing matrix information to achieve increased efficiency. 3752 3753 Options Database Keys: 3754 + -mat_no_inode - Do not use inodes 3755 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3756 3757 Level: intermediate 3758 3759 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3760 3761 @*/ 3762 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3763 { 3764 PetscErrorCode ierr; 3765 3766 PetscFunctionBegin; 3767 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3768 PetscValidType(B,1); 3769 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3770 PetscFunctionReturn(0); 3771 } 3772 3773 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3774 { 3775 Mat_SeqAIJ *b; 3776 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3777 PetscErrorCode ierr; 3778 PetscInt i; 3779 3780 PetscFunctionBegin; 3781 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3782 if (nz == MAT_SKIP_ALLOCATION) { 3783 skipallocation = PETSC_TRUE; 3784 nz = 0; 3785 } 3786 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3787 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3788 3789 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3790 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3791 if (nnz) { 3792 for (i=0; i<B->rmap->n; i++) { 3793 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]); 3794 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %D value %d rowlength %D",i,nnz[i],B->cmap->n); 3795 } 3796 } 3797 3798 B->preallocated = PETSC_TRUE; 3799 3800 b = (Mat_SeqAIJ*)B->data; 3801 3802 if (!skipallocation) { 3803 if (!b->imax) { 3804 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3805 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3806 } 3807 if (!b->ipre) { 3808 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3809 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3810 } 3811 if (!nnz) { 3812 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3813 else if (nz < 0) nz = 1; 3814 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3815 nz = nz*B->rmap->n; 3816 } else { 3817 nz = 0; 3818 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3819 } 3820 /* b->ilen will count nonzeros in each row so far. */ 3821 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3822 3823 /* allocate the matrix space */ 3824 /* FIXME: should B's old memory be unlogged? */ 3825 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3826 if (B->structure_only) { 3827 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3828 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3829 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3830 } else { 3831 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3832 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3833 } 3834 b->i[0] = 0; 3835 for (i=1; i<B->rmap->n+1; i++) { 3836 b->i[i] = b->i[i-1] + b->imax[i-1]; 3837 } 3838 if (B->structure_only) { 3839 b->singlemalloc = PETSC_FALSE; 3840 b->free_a = PETSC_FALSE; 3841 } else { 3842 b->singlemalloc = PETSC_TRUE; 3843 b->free_a = PETSC_TRUE; 3844 } 3845 b->free_ij = PETSC_TRUE; 3846 } else { 3847 b->free_a = PETSC_FALSE; 3848 b->free_ij = PETSC_FALSE; 3849 } 3850 3851 if (b->ipre && nnz != b->ipre && b->imax) { 3852 /* reserve user-requested sparsity */ 3853 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3854 } 3855 3856 3857 b->nz = 0; 3858 b->maxnz = nz; 3859 B->info.nz_unneeded = (double)b->maxnz; 3860 if (realalloc) { 3861 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3862 } 3863 B->was_assembled = PETSC_FALSE; 3864 B->assembled = PETSC_FALSE; 3865 PetscFunctionReturn(0); 3866 } 3867 3868 3869 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3870 { 3871 Mat_SeqAIJ *a; 3872 PetscInt i; 3873 PetscErrorCode ierr; 3874 3875 PetscFunctionBegin; 3876 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3877 a = (Mat_SeqAIJ*)A->data; 3878 /* if no saved info, we error out */ 3879 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3880 3881 if (!a->i || !a->j || !a->a || !a->imax || !a->ilen) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation \n"); 3882 3883 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3884 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3885 a->i[0] = 0; 3886 for (i=1; i<A->rmap->n+1; i++) { 3887 a->i[i] = a->i[i-1] + a->imax[i-1]; 3888 } 3889 A->preallocated = PETSC_TRUE; 3890 a->nz = 0; 3891 a->maxnz = a->i[A->rmap->n]; 3892 A->info.nz_unneeded = (double)a->maxnz; 3893 A->was_assembled = PETSC_FALSE; 3894 A->assembled = PETSC_FALSE; 3895 PetscFunctionReturn(0); 3896 } 3897 3898 /*@ 3899 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3900 3901 Input Parameters: 3902 + B - the matrix 3903 . i - the indices into j for the start of each row (starts with zero) 3904 . j - the column indices for each row (starts with zero) these must be sorted for each row 3905 - v - optional values in the matrix 3906 3907 Level: developer 3908 3909 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3910 3911 .keywords: matrix, aij, compressed row, sparse, sequential 3912 3913 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ 3914 @*/ 3915 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3916 { 3917 PetscErrorCode ierr; 3918 3919 PetscFunctionBegin; 3920 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3921 PetscValidType(B,1); 3922 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3923 PetscFunctionReturn(0); 3924 } 3925 3926 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3927 { 3928 PetscInt i; 3929 PetscInt m,n; 3930 PetscInt nz; 3931 PetscInt *nnz, nz_max = 0; 3932 PetscScalar *values; 3933 PetscErrorCode ierr; 3934 3935 PetscFunctionBegin; 3936 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3937 3938 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3939 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3940 3941 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3942 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3943 for (i = 0; i < m; i++) { 3944 nz = Ii[i+1]- Ii[i]; 3945 nz_max = PetscMax(nz_max, nz); 3946 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3947 nnz[i] = nz; 3948 } 3949 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3950 ierr = PetscFree(nnz);CHKERRQ(ierr); 3951 3952 if (v) { 3953 values = (PetscScalar*) v; 3954 } else { 3955 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3956 } 3957 3958 for (i = 0; i < m; i++) { 3959 nz = Ii[i+1] - Ii[i]; 3960 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3961 } 3962 3963 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3964 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3965 3966 if (!v) { 3967 ierr = PetscFree(values);CHKERRQ(ierr); 3968 } 3969 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3970 PetscFunctionReturn(0); 3971 } 3972 3973 #include <../src/mat/impls/dense/seq/dense.h> 3974 #include <petsc/private/kernels/petscaxpy.h> 3975 3976 /* 3977 Computes (B'*A')' since computing B*A directly is untenable 3978 3979 n p p 3980 ( ) ( ) ( ) 3981 m ( A ) * n ( B ) = m ( C ) 3982 ( ) ( ) ( ) 3983 3984 */ 3985 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3986 { 3987 PetscErrorCode ierr; 3988 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3989 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3990 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3991 PetscInt i,n,m,q,p; 3992 const PetscInt *ii,*idx; 3993 const PetscScalar *b,*a,*a_q; 3994 PetscScalar *c,*c_q; 3995 3996 PetscFunctionBegin; 3997 m = A->rmap->n; 3998 n = A->cmap->n; 3999 p = B->cmap->n; 4000 a = sub_a->v; 4001 b = sub_b->a; 4002 c = sub_c->v; 4003 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 4004 4005 ii = sub_b->i; 4006 idx = sub_b->j; 4007 for (i=0; i<n; i++) { 4008 q = ii[i+1] - ii[i]; 4009 while (q-->0) { 4010 c_q = c + m*(*idx); 4011 a_q = a + m*i; 4012 PetscKernelAXPY(c_q,*b,a_q,m); 4013 idx++; 4014 b++; 4015 } 4016 } 4017 PetscFunctionReturn(0); 4018 } 4019 4020 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 4021 { 4022 PetscErrorCode ierr; 4023 PetscInt m=A->rmap->n,n=B->cmap->n; 4024 Mat Cmat; 4025 4026 PetscFunctionBegin; 4027 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %D != B->rmap->n %D\n",A->cmap->n,B->rmap->n); 4028 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 4029 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 4030 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 4031 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 4032 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 4033 4034 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4035 4036 *C = Cmat; 4037 PetscFunctionReturn(0); 4038 } 4039 4040 /* ----------------------------------------------------------------*/ 4041 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 4042 { 4043 PetscErrorCode ierr; 4044 4045 PetscFunctionBegin; 4046 if (scall == MAT_INITIAL_MATRIX) { 4047 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4048 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 4049 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4050 } 4051 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4052 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 4053 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4054 PetscFunctionReturn(0); 4055 } 4056 4057 4058 /*MC 4059 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4060 based on compressed sparse row format. 4061 4062 Options Database Keys: 4063 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4064 4065 Level: beginner 4066 4067 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 4068 M*/ 4069 4070 /*MC 4071 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4072 4073 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4074 and MATMPIAIJ otherwise. As a result, for single process communicators, 4075 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 4076 for communicators controlling multiple processes. It is recommended that you call both of 4077 the above preallocation routines for simplicity. 4078 4079 Options Database Keys: 4080 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4081 4082 Developer Notes: 4083 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4084 enough exist. 4085 4086 Level: beginner 4087 4088 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 4089 M*/ 4090 4091 /*MC 4092 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4093 4094 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4095 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4096 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4097 for communicators controlling multiple processes. It is recommended that you call both of 4098 the above preallocation routines for simplicity. 4099 4100 Options Database Keys: 4101 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4102 4103 Level: beginner 4104 4105 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4106 M*/ 4107 4108 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4109 #if defined(PETSC_HAVE_ELEMENTAL) 4110 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4111 #endif 4112 #if defined(PETSC_HAVE_HYPRE) 4113 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4114 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 4115 #endif 4116 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 4117 4118 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4119 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4120 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 4121 4122 /*@C 4123 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 4124 4125 Not Collective 4126 4127 Input Parameter: 4128 . mat - a MATSEQAIJ matrix 4129 4130 Output Parameter: 4131 . array - pointer to the data 4132 4133 Level: intermediate 4134 4135 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4136 @*/ 4137 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4138 { 4139 PetscErrorCode ierr; 4140 4141 PetscFunctionBegin; 4142 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4143 PetscFunctionReturn(0); 4144 } 4145 4146 /*@C 4147 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4148 4149 Not Collective 4150 4151 Input Parameter: 4152 . mat - a MATSEQAIJ matrix 4153 4154 Output Parameter: 4155 . nz - the maximum number of nonzeros in any row 4156 4157 Level: intermediate 4158 4159 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4160 @*/ 4161 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4162 { 4163 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4164 4165 PetscFunctionBegin; 4166 *nz = aij->rmax; 4167 PetscFunctionReturn(0); 4168 } 4169 4170 /*@C 4171 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4172 4173 Not Collective 4174 4175 Input Parameters: 4176 . mat - a MATSEQAIJ matrix 4177 . array - pointer to the data 4178 4179 Level: intermediate 4180 4181 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4182 @*/ 4183 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4184 { 4185 PetscErrorCode ierr; 4186 4187 PetscFunctionBegin; 4188 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4189 PetscFunctionReturn(0); 4190 } 4191 4192 #if defined(PETSC_HAVE_CUDA) 4193 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat); 4194 #endif 4195 4196 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4197 { 4198 Mat_SeqAIJ *b; 4199 PetscErrorCode ierr; 4200 PetscMPIInt size; 4201 4202 PetscFunctionBegin; 4203 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4204 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4205 4206 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4207 4208 B->data = (void*)b; 4209 4210 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4211 4212 b->row = 0; 4213 b->col = 0; 4214 b->icol = 0; 4215 b->reallocs = 0; 4216 b->ignorezeroentries = PETSC_FALSE; 4217 b->roworiented = PETSC_TRUE; 4218 b->nonew = 0; 4219 b->diag = 0; 4220 b->solve_work = 0; 4221 B->spptr = 0; 4222 b->saved_values = 0; 4223 b->idiag = 0; 4224 b->mdiag = 0; 4225 b->ssor_work = 0; 4226 b->omega = 1.0; 4227 b->fshift = 0.0; 4228 b->idiagvalid = PETSC_FALSE; 4229 b->ibdiagvalid = PETSC_FALSE; 4230 b->keepnonzeropattern = PETSC_FALSE; 4231 4232 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4233 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4234 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4235 4236 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4237 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4238 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4239 #endif 4240 4241 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4242 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4243 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4244 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4245 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4246 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4247 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4248 #if defined(PETSC_HAVE_MKL_SPARSE) 4249 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4250 #endif 4251 #if defined(PETSC_HAVE_CUDA) 4252 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);CHKERRQ(ierr); 4253 #endif 4254 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4255 #if defined(PETSC_HAVE_ELEMENTAL) 4256 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4257 #endif 4258 #if defined(PETSC_HAVE_HYPRE) 4259 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4260 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4261 #endif 4262 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4263 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4264 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4265 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4266 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4267 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4268 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4269 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4270 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4271 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4272 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4273 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4274 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaij_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 4275 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4276 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4277 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4278 PetscFunctionReturn(0); 4279 } 4280 4281 /* 4282 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4283 */ 4284 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4285 { 4286 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4287 PetscErrorCode ierr; 4288 PetscInt i,m = A->rmap->n; 4289 4290 PetscFunctionBegin; 4291 c = (Mat_SeqAIJ*)C->data; 4292 4293 C->factortype = A->factortype; 4294 c->row = 0; 4295 c->col = 0; 4296 c->icol = 0; 4297 c->reallocs = 0; 4298 4299 C->assembled = PETSC_TRUE; 4300 4301 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4302 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4303 4304 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4305 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4306 for (i=0; i<m; i++) { 4307 c->imax[i] = a->imax[i]; 4308 c->ilen[i] = a->ilen[i]; 4309 } 4310 4311 /* allocate the matrix space */ 4312 if (mallocmatspace) { 4313 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4314 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4315 4316 c->singlemalloc = PETSC_TRUE; 4317 4318 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4319 if (m > 0) { 4320 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4321 if (cpvalues == MAT_COPY_VALUES) { 4322 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4323 } else { 4324 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4325 } 4326 } 4327 } 4328 4329 c->ignorezeroentries = a->ignorezeroentries; 4330 c->roworiented = a->roworiented; 4331 c->nonew = a->nonew; 4332 if (a->diag) { 4333 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4334 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4335 for (i=0; i<m; i++) { 4336 c->diag[i] = a->diag[i]; 4337 } 4338 } else c->diag = 0; 4339 4340 c->solve_work = 0; 4341 c->saved_values = 0; 4342 c->idiag = 0; 4343 c->ssor_work = 0; 4344 c->keepnonzeropattern = a->keepnonzeropattern; 4345 c->free_a = PETSC_TRUE; 4346 c->free_ij = PETSC_TRUE; 4347 4348 c->rmax = a->rmax; 4349 c->nz = a->nz; 4350 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4351 C->preallocated = PETSC_TRUE; 4352 4353 c->compressedrow.use = a->compressedrow.use; 4354 c->compressedrow.nrows = a->compressedrow.nrows; 4355 if (a->compressedrow.use) { 4356 i = a->compressedrow.nrows; 4357 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4358 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4359 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4360 } else { 4361 c->compressedrow.use = PETSC_FALSE; 4362 c->compressedrow.i = NULL; 4363 c->compressedrow.rindex = NULL; 4364 } 4365 c->nonzerorowcnt = a->nonzerorowcnt; 4366 C->nonzerostate = A->nonzerostate; 4367 4368 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4369 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4370 PetscFunctionReturn(0); 4371 } 4372 4373 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4374 { 4375 PetscErrorCode ierr; 4376 4377 PetscFunctionBegin; 4378 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4379 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4380 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4381 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4382 } 4383 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4384 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4385 PetscFunctionReturn(0); 4386 } 4387 4388 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4389 { 4390 PetscBool isbinary, ishdf5; 4391 PetscErrorCode ierr; 4392 4393 PetscFunctionBegin; 4394 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4395 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4396 /* force binary viewer to load .info file if it has not yet done so */ 4397 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4398 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 4399 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 4400 if (isbinary) { 4401 ierr = MatLoad_SeqAIJ_Binary(newMat,viewer);CHKERRQ(ierr); 4402 } else if (ishdf5) { 4403 #if defined(PETSC_HAVE_HDF5) 4404 ierr = MatLoad_AIJ_HDF5(newMat,viewer);CHKERRQ(ierr); 4405 #else 4406 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4407 #endif 4408 } else { 4409 SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 4410 } 4411 PetscFunctionReturn(0); 4412 } 4413 4414 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat newMat, PetscViewer viewer) 4415 { 4416 Mat_SeqAIJ *a; 4417 PetscErrorCode ierr; 4418 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4419 int fd; 4420 PetscMPIInt size; 4421 MPI_Comm comm; 4422 PetscInt bs = newMat->rmap->bs; 4423 4424 PetscFunctionBegin; 4425 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4426 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4427 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4428 4429 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4430 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4431 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4432 if (bs < 0) bs = 1; 4433 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4434 4435 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4436 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4437 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4438 M = header[1]; N = header[2]; nz = header[3]; 4439 4440 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4441 4442 /* read in row lengths */ 4443 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4444 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4445 4446 /* check if sum of rowlengths is same as nz */ 4447 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4448 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %dD, sum-row-lengths = %D\n",nz,sum); 4449 4450 /* set global size if not set already*/ 4451 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4452 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4453 } else { 4454 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4455 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4456 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4457 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4458 } 4459 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); 4460 } 4461 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4462 a = (Mat_SeqAIJ*)newMat->data; 4463 4464 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4465 4466 /* read in nonzero values */ 4467 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4468 4469 /* set matrix "i" values */ 4470 a->i[0] = 0; 4471 for (i=1; i<= M; i++) { 4472 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4473 a->ilen[i-1] = rowlengths[i-1]; 4474 } 4475 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4476 4477 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4478 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4479 PetscFunctionReturn(0); 4480 } 4481 4482 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4483 { 4484 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4485 PetscErrorCode ierr; 4486 #if defined(PETSC_USE_COMPLEX) 4487 PetscInt k; 4488 #endif 4489 4490 PetscFunctionBegin; 4491 /* If the matrix dimensions are not equal,or no of nonzeros */ 4492 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4493 *flg = PETSC_FALSE; 4494 PetscFunctionReturn(0); 4495 } 4496 4497 /* if the a->i are the same */ 4498 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4499 if (!*flg) PetscFunctionReturn(0); 4500 4501 /* if a->j are the same */ 4502 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4503 if (!*flg) PetscFunctionReturn(0); 4504 4505 /* if a->a are the same */ 4506 #if defined(PETSC_USE_COMPLEX) 4507 for (k=0; k<a->nz; k++) { 4508 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4509 *flg = PETSC_FALSE; 4510 PetscFunctionReturn(0); 4511 } 4512 } 4513 #else 4514 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4515 #endif 4516 PetscFunctionReturn(0); 4517 } 4518 4519 /*@ 4520 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4521 provided by the user. 4522 4523 Collective on MPI_Comm 4524 4525 Input Parameters: 4526 + comm - must be an MPI communicator of size 1 4527 . m - number of rows 4528 . n - number of columns 4529 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 4530 . j - column indices 4531 - a - matrix values 4532 4533 Output Parameter: 4534 . mat - the matrix 4535 4536 Level: intermediate 4537 4538 Notes: 4539 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4540 once the matrix is destroyed and not before 4541 4542 You cannot set new nonzero locations into this matrix, that will generate an error. 4543 4544 The i and j indices are 0 based 4545 4546 The format which is used for the sparse matrix input, is equivalent to a 4547 row-major ordering.. i.e for the following matrix, the input data expected is 4548 as shown 4549 4550 $ 1 0 0 4551 $ 2 0 3 4552 $ 4 5 6 4553 $ 4554 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4555 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4556 $ v = {1,2,3,4,5,6} [size = 6] 4557 4558 4559 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4560 4561 @*/ 4562 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4563 { 4564 PetscErrorCode ierr; 4565 PetscInt ii; 4566 Mat_SeqAIJ *aij; 4567 #if defined(PETSC_USE_DEBUG) 4568 PetscInt jj; 4569 #endif 4570 4571 PetscFunctionBegin; 4572 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4573 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4574 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4575 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4576 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4577 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4578 aij = (Mat_SeqAIJ*)(*mat)->data; 4579 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4580 4581 aij->i = i; 4582 aij->j = j; 4583 aij->a = a; 4584 aij->singlemalloc = PETSC_FALSE; 4585 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4586 aij->free_a = PETSC_FALSE; 4587 aij->free_ij = PETSC_FALSE; 4588 4589 for (ii=0; ii<m; ii++) { 4590 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4591 #if defined(PETSC_USE_DEBUG) 4592 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]); 4593 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4594 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual column %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 4595 if (j[jj] == j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual column %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 4596 } 4597 #endif 4598 } 4599 #if defined(PETSC_USE_DEBUG) 4600 for (ii=0; ii<aij->i[m]; ii++) { 4601 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4602 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]); 4603 } 4604 #endif 4605 4606 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4607 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4608 PetscFunctionReturn(0); 4609 } 4610 /*@C 4611 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4612 provided by the user. 4613 4614 Collective on MPI_Comm 4615 4616 Input Parameters: 4617 + comm - must be an MPI communicator of size 1 4618 . m - number of rows 4619 . n - number of columns 4620 . i - row indices 4621 . j - column indices 4622 . a - matrix values 4623 . nz - number of nonzeros 4624 - idx - 0 or 1 based 4625 4626 Output Parameter: 4627 . mat - the matrix 4628 4629 Level: intermediate 4630 4631 Notes: 4632 The i and j indices are 0 based 4633 4634 The format which is used for the sparse matrix input, is equivalent to a 4635 row-major ordering.. i.e for the following matrix, the input data expected is 4636 as shown: 4637 4638 1 0 0 4639 2 0 3 4640 4 5 6 4641 4642 i = {0,1,1,2,2,2} 4643 j = {0,0,2,0,1,2} 4644 v = {1,2,3,4,5,6} 4645 4646 4647 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4648 4649 @*/ 4650 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4651 { 4652 PetscErrorCode ierr; 4653 PetscInt ii, *nnz, one = 1,row,col; 4654 4655 4656 PetscFunctionBegin; 4657 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4658 for (ii = 0; ii < nz; ii++) { 4659 nnz[i[ii] - !!idx] += 1; 4660 } 4661 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4662 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4663 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4664 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4665 for (ii = 0; ii < nz; ii++) { 4666 if (idx) { 4667 row = i[ii] - 1; 4668 col = j[ii] - 1; 4669 } else { 4670 row = i[ii]; 4671 col = j[ii]; 4672 } 4673 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4674 } 4675 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4676 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4677 ierr = PetscFree(nnz);CHKERRQ(ierr); 4678 PetscFunctionReturn(0); 4679 } 4680 4681 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4682 { 4683 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4684 PetscErrorCode ierr; 4685 4686 PetscFunctionBegin; 4687 a->idiagvalid = PETSC_FALSE; 4688 a->ibdiagvalid = PETSC_FALSE; 4689 4690 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4691 PetscFunctionReturn(0); 4692 } 4693 4694 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4695 { 4696 PetscErrorCode ierr; 4697 PetscMPIInt size; 4698 4699 PetscFunctionBegin; 4700 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4701 if (size == 1) { 4702 if (scall == MAT_INITIAL_MATRIX) { 4703 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4704 } else { 4705 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4706 } 4707 } else { 4708 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4709 } 4710 PetscFunctionReturn(0); 4711 } 4712 4713 /* 4714 Permute A into C's *local* index space using rowemb,colemb. 4715 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4716 of [0,m), colemb is in [0,n). 4717 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4718 */ 4719 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4720 { 4721 /* If making this function public, change the error returned in this function away from _PLIB. */ 4722 PetscErrorCode ierr; 4723 Mat_SeqAIJ *Baij; 4724 PetscBool seqaij; 4725 PetscInt m,n,*nz,i,j,count; 4726 PetscScalar v; 4727 const PetscInt *rowindices,*colindices; 4728 4729 PetscFunctionBegin; 4730 if (!B) PetscFunctionReturn(0); 4731 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4732 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4733 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4734 if (rowemb) { 4735 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4736 if (m != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %D is incompatible with matrix row size %D",m,B->rmap->n); 4737 } else { 4738 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4739 } 4740 if (colemb) { 4741 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4742 if (n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %D is incompatible with input matrix col size %D",n,B->cmap->n); 4743 } else { 4744 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4745 } 4746 4747 Baij = (Mat_SeqAIJ*)(B->data); 4748 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4749 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4750 for (i=0; i<B->rmap->n; i++) { 4751 nz[i] = Baij->i[i+1] - Baij->i[i]; 4752 } 4753 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4754 ierr = PetscFree(nz);CHKERRQ(ierr); 4755 } 4756 if (pattern == SUBSET_NONZERO_PATTERN) { 4757 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4758 } 4759 count = 0; 4760 rowindices = NULL; 4761 colindices = NULL; 4762 if (rowemb) { 4763 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4764 } 4765 if (colemb) { 4766 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4767 } 4768 for (i=0; i<B->rmap->n; i++) { 4769 PetscInt row; 4770 row = i; 4771 if (rowindices) row = rowindices[i]; 4772 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4773 PetscInt col; 4774 col = Baij->j[count]; 4775 if (colindices) col = colindices[col]; 4776 v = Baij->a[count]; 4777 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4778 ++count; 4779 } 4780 } 4781 /* FIXME: set C's nonzerostate correctly. */ 4782 /* Assembly for C is necessary. */ 4783 C->preallocated = PETSC_TRUE; 4784 C->assembled = PETSC_TRUE; 4785 C->was_assembled = PETSC_FALSE; 4786 PetscFunctionReturn(0); 4787 } 4788 4789 PetscFunctionList MatSeqAIJList = NULL; 4790 4791 /*@C 4792 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4793 4794 Collective on Mat 4795 4796 Input Parameters: 4797 + mat - the matrix object 4798 - matype - matrix type 4799 4800 Options Database Key: 4801 . -mat_seqai_type <method> - for example seqaijcrl 4802 4803 4804 Level: intermediate 4805 4806 .keywords: Mat, MatType, set, method 4807 4808 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4809 @*/ 4810 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4811 { 4812 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4813 PetscBool sametype; 4814 4815 PetscFunctionBegin; 4816 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4817 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4818 if (sametype) PetscFunctionReturn(0); 4819 4820 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4821 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4822 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4823 PetscFunctionReturn(0); 4824 } 4825 4826 4827 /*@C 4828 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4829 4830 Not Collective 4831 4832 Input Parameters: 4833 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4834 - function - routine to convert to subtype 4835 4836 Notes: 4837 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4838 4839 4840 Then, your matrix can be chosen with the procedural interface at runtime via the option 4841 $ -mat_seqaij_type my_mat 4842 4843 Level: advanced 4844 4845 .keywords: Mat, register 4846 4847 .seealso: MatSeqAIJRegisterAll() 4848 4849 4850 Level: advanced 4851 @*/ 4852 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4853 { 4854 PetscErrorCode ierr; 4855 4856 PetscFunctionBegin; 4857 ierr = MatInitializePackage();CHKERRQ(ierr); 4858 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4859 PetscFunctionReturn(0); 4860 } 4861 4862 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4863 4864 /*@C 4865 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4866 4867 Not Collective 4868 4869 Level: advanced 4870 4871 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4872 4873 .keywords: KSP, register, all 4874 4875 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4876 @*/ 4877 PetscErrorCode MatSeqAIJRegisterAll(void) 4878 { 4879 PetscErrorCode ierr; 4880 4881 PetscFunctionBegin; 4882 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4883 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4884 4885 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4886 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4887 ierr = MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4888 #if defined(PETSC_HAVE_MKL_SPARSE) 4889 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4890 #endif 4891 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4892 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4893 #endif 4894 PetscFunctionReturn(0); 4895 } 4896 4897 /* 4898 Special version for direct calls from Fortran 4899 */ 4900 #include <petsc/private/fortranimpl.h> 4901 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4902 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4903 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4904 #define matsetvaluesseqaij_ matsetvaluesseqaij 4905 #endif 4906 4907 /* Change these macros so can be used in void function */ 4908 #undef CHKERRQ 4909 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4910 #undef SETERRQ2 4911 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4912 #undef SETERRQ3 4913 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4914 4915 PETSC_EXTERN void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4916 { 4917 Mat A = *AA; 4918 PetscInt m = *mm, n = *nn; 4919 InsertMode is = *isis; 4920 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4921 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4922 PetscInt *imax,*ai,*ailen; 4923 PetscErrorCode ierr; 4924 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4925 MatScalar *ap,value,*aa; 4926 PetscBool ignorezeroentries = a->ignorezeroentries; 4927 PetscBool roworiented = a->roworiented; 4928 4929 PetscFunctionBegin; 4930 MatCheckPreallocated(A,1); 4931 imax = a->imax; 4932 ai = a->i; 4933 ailen = a->ilen; 4934 aj = a->j; 4935 aa = a->a; 4936 4937 for (k=0; k<m; k++) { /* loop over added rows */ 4938 row = im[k]; 4939 if (row < 0) continue; 4940 #if defined(PETSC_USE_DEBUG) 4941 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4942 #endif 4943 rp = aj + ai[row]; ap = aa + ai[row]; 4944 rmax = imax[row]; nrow = ailen[row]; 4945 low = 0; 4946 high = nrow; 4947 for (l=0; l<n; l++) { /* loop over added columns */ 4948 if (in[l] < 0) continue; 4949 #if defined(PETSC_USE_DEBUG) 4950 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4951 #endif 4952 col = in[l]; 4953 if (roworiented) value = v[l + k*n]; 4954 else value = v[k + l*m]; 4955 4956 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4957 4958 if (col <= lastcol) low = 0; 4959 else high = nrow; 4960 lastcol = col; 4961 while (high-low > 5) { 4962 t = (low+high)/2; 4963 if (rp[t] > col) high = t; 4964 else low = t; 4965 } 4966 for (i=low; i<high; i++) { 4967 if (rp[i] > col) break; 4968 if (rp[i] == col) { 4969 if (is == ADD_VALUES) ap[i] += value; 4970 else ap[i] = value; 4971 goto noinsert; 4972 } 4973 } 4974 if (value == 0.0 && ignorezeroentries) goto noinsert; 4975 if (nonew == 1) goto noinsert; 4976 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4977 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4978 N = nrow++ - 1; a->nz++; high++; 4979 /* shift up all the later entries in this row */ 4980 for (ii=N; ii>=i; ii--) { 4981 rp[ii+1] = rp[ii]; 4982 ap[ii+1] = ap[ii]; 4983 } 4984 rp[i] = col; 4985 ap[i] = value; 4986 A->nonzerostate++; 4987 noinsert:; 4988 low = i + 1; 4989 } 4990 ailen[row] = nrow; 4991 } 4992 PetscFunctionReturnVoid(); 4993 } 4994