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 a->imax[i] = PetscMin(a->imax[i],A->cmap->n); 1563 } 1564 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax);CHKERRQ(ierr); 1565 1566 /* copy old values into new matrix data structure */ 1567 for (i=0; i<A->rmap->n; i++) { 1568 ierr = MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES);CHKERRQ(ierr); 1569 if (i < A->cmap->n) { 1570 ierr = MatSetValue(A,i,i,v,ADD_VALUES);CHKERRQ(ierr); 1571 } 1572 } 1573 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1574 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1575 if (singlemalloc) { 1576 ierr = PetscFree3(olda,oldj,oldi);CHKERRQ(ierr); 1577 } else { 1578 if (free_a) {ierr = PetscFree(olda);CHKERRQ(ierr);} 1579 if (free_ij) {ierr = PetscFree(oldj);CHKERRQ(ierr);} 1580 if (free_ij) {ierr = PetscFree(oldi);CHKERRQ(ierr);} 1581 } 1582 } 1583 ierr = PetscFree(mdiag);CHKERRQ(ierr); 1584 a->diagonaldense = PETSC_TRUE; 1585 PetscFunctionReturn(0); 1586 } 1587 1588 /* 1589 Checks for missing diagonals 1590 */ 1591 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1592 { 1593 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1594 PetscInt *diag,*ii = a->i,i; 1595 PetscErrorCode ierr; 1596 1597 PetscFunctionBegin; 1598 *missing = PETSC_FALSE; 1599 if (A->rmap->n > 0 && !ii) { 1600 *missing = PETSC_TRUE; 1601 if (d) *d = 0; 1602 ierr = PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");CHKERRQ(ierr); 1603 } else { 1604 diag = a->diag; 1605 for (i=0; i<A->rmap->n; i++) { 1606 if (diag[i] >= ii[i+1]) { 1607 *missing = PETSC_TRUE; 1608 if (d) *d = i; 1609 ierr = PetscInfo1(A,"Matrix is missing diagonal number %D\n",i);CHKERRQ(ierr); 1610 break; 1611 } 1612 } 1613 } 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #include <petscblaslapack.h> 1618 #include <petsc/private/kernels/blockinvert.h> 1619 1620 /* 1621 Note that values is allocated externally by the PC and then passed into this routine 1622 */ 1623 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1624 { 1625 PetscErrorCode ierr; 1626 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1627 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1628 const PetscReal shift = 0.0; 1629 PetscInt ipvt[5]; 1630 PetscScalar work[25],*v_work; 1631 1632 PetscFunctionBegin; 1633 allowzeropivot = PetscNot(A->erroriffailure); 1634 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1635 if (ncnt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %D doesn't match number matrix rows %D",ncnt,n); 1636 for (i=0; i<nblocks; i++) { 1637 bsizemax = PetscMax(bsizemax,bsizes[i]); 1638 } 1639 ierr = PetscMalloc1(bsizemax,&indx);CHKERRQ(ierr); 1640 if (bsizemax > 7) { 1641 ierr = PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots);CHKERRQ(ierr); 1642 } 1643 ncnt = 0; 1644 for (i=0; i<nblocks; i++) { 1645 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1646 ierr = MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag);CHKERRQ(ierr); 1647 switch (bsizes[i]) { 1648 case 1: 1649 *diag = 1.0/(*diag); 1650 break; 1651 case 2: 1652 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1653 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1654 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 1655 break; 1656 case 3: 1657 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1658 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1659 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 1660 break; 1661 case 4: 1662 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1663 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1664 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 1665 break; 1666 case 5: 1667 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1668 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1669 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 1670 break; 1671 case 6: 1672 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1673 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1674 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 1675 break; 1676 case 7: 1677 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1678 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1679 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 1680 break; 1681 default: 1682 ierr = PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1683 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1684 ierr = PetscKernel_A_gets_transpose_A_N(diag,bsizes[i]);CHKERRQ(ierr); 1685 } 1686 ncnt += bsizes[i]; 1687 diag += bsizes[i]*bsizes[i]; 1688 } 1689 if (bsizemax > 7) { 1690 ierr = PetscFree2(v_work,v_pivots);CHKERRQ(ierr); 1691 } 1692 ierr = PetscFree(indx);CHKERRQ(ierr); 1693 PetscFunctionReturn(0); 1694 } 1695 1696 /* 1697 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1698 */ 1699 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1700 { 1701 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1702 PetscErrorCode ierr; 1703 PetscInt i,*diag,m = A->rmap->n; 1704 MatScalar *v = a->a; 1705 PetscScalar *idiag,*mdiag; 1706 1707 PetscFunctionBegin; 1708 if (a->idiagvalid) PetscFunctionReturn(0); 1709 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1710 diag = a->diag; 1711 if (!a->idiag) { 1712 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1713 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1714 v = a->a; 1715 } 1716 mdiag = a->mdiag; 1717 idiag = a->idiag; 1718 1719 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1720 for (i=0; i<m; i++) { 1721 mdiag[i] = v[diag[i]]; 1722 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1723 if (PetscRealPart(fshift)) { 1724 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1725 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1726 A->factorerror_zeropivot_value = 0.0; 1727 A->factorerror_zeropivot_row = i; 1728 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1729 } 1730 idiag[i] = 1.0/v[diag[i]]; 1731 } 1732 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1733 } else { 1734 for (i=0; i<m; i++) { 1735 mdiag[i] = v[diag[i]]; 1736 idiag[i] = omega/(fshift + v[diag[i]]); 1737 } 1738 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1739 } 1740 a->idiagvalid = PETSC_TRUE; 1741 PetscFunctionReturn(0); 1742 } 1743 1744 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1745 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1746 { 1747 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1748 PetscScalar *x,d,sum,*t,scale; 1749 const MatScalar *v,*idiag=0,*mdiag; 1750 const PetscScalar *b, *bs,*xb, *ts; 1751 PetscErrorCode ierr; 1752 PetscInt n,m = A->rmap->n,i; 1753 const PetscInt *idx,*diag; 1754 1755 PetscFunctionBegin; 1756 its = its*lits; 1757 1758 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1759 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1760 a->fshift = fshift; 1761 a->omega = omega; 1762 1763 diag = a->diag; 1764 t = a->ssor_work; 1765 idiag = a->idiag; 1766 mdiag = a->mdiag; 1767 1768 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1769 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1770 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1771 if (flag == SOR_APPLY_UPPER) { 1772 /* apply (U + D/omega) to the vector */ 1773 bs = b; 1774 for (i=0; i<m; i++) { 1775 d = fshift + mdiag[i]; 1776 n = a->i[i+1] - diag[i] - 1; 1777 idx = a->j + diag[i] + 1; 1778 v = a->a + diag[i] + 1; 1779 sum = b[i]*d/omega; 1780 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1781 x[i] = sum; 1782 } 1783 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1784 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1785 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1786 PetscFunctionReturn(0); 1787 } 1788 1789 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1790 else if (flag & SOR_EISENSTAT) { 1791 /* Let A = L + U + D; where L is lower trianglar, 1792 U is upper triangular, E = D/omega; This routine applies 1793 1794 (L + E)^{-1} A (U + E)^{-1} 1795 1796 to a vector efficiently using Eisenstat's trick. 1797 */ 1798 scale = (2.0/omega) - 1.0; 1799 1800 /* x = (E + U)^{-1} b */ 1801 for (i=m-1; i>=0; i--) { 1802 n = a->i[i+1] - diag[i] - 1; 1803 idx = a->j + diag[i] + 1; 1804 v = a->a + diag[i] + 1; 1805 sum = b[i]; 1806 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1807 x[i] = sum*idiag[i]; 1808 } 1809 1810 /* t = b - (2*E - D)x */ 1811 v = a->a; 1812 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1813 1814 /* t = (E + L)^{-1}t */ 1815 ts = t; 1816 diag = a->diag; 1817 for (i=0; i<m; i++) { 1818 n = diag[i] - a->i[i]; 1819 idx = a->j + a->i[i]; 1820 v = a->a + a->i[i]; 1821 sum = t[i]; 1822 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1823 t[i] = sum*idiag[i]; 1824 /* x = x + t */ 1825 x[i] += t[i]; 1826 } 1827 1828 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1829 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1830 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1831 PetscFunctionReturn(0); 1832 } 1833 if (flag & SOR_ZERO_INITIAL_GUESS) { 1834 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1835 for (i=0; i<m; i++) { 1836 n = diag[i] - a->i[i]; 1837 idx = a->j + a->i[i]; 1838 v = a->a + a->i[i]; 1839 sum = b[i]; 1840 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1841 t[i] = sum; 1842 x[i] = sum*idiag[i]; 1843 } 1844 xb = t; 1845 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1846 } else xb = b; 1847 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1848 for (i=m-1; i>=0; i--) { 1849 n = a->i[i+1] - diag[i] - 1; 1850 idx = a->j + diag[i] + 1; 1851 v = a->a + diag[i] + 1; 1852 sum = xb[i]; 1853 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1854 if (xb == b) { 1855 x[i] = sum*idiag[i]; 1856 } else { 1857 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1858 } 1859 } 1860 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1861 } 1862 its--; 1863 } 1864 while (its--) { 1865 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1866 for (i=0; i<m; i++) { 1867 /* lower */ 1868 n = diag[i] - a->i[i]; 1869 idx = a->j + a->i[i]; 1870 v = a->a + a->i[i]; 1871 sum = b[i]; 1872 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1873 t[i] = sum; /* save application of the lower-triangular part */ 1874 /* upper */ 1875 n = a->i[i+1] - diag[i] - 1; 1876 idx = a->j + diag[i] + 1; 1877 v = a->a + diag[i] + 1; 1878 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1879 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1880 } 1881 xb = t; 1882 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1883 } else xb = b; 1884 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1885 for (i=m-1; i>=0; i--) { 1886 sum = xb[i]; 1887 if (xb == b) { 1888 /* whole matrix (no checkpointing available) */ 1889 n = a->i[i+1] - a->i[i]; 1890 idx = a->j + a->i[i]; 1891 v = a->a + a->i[i]; 1892 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1893 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1894 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1895 n = a->i[i+1] - diag[i] - 1; 1896 idx = a->j + diag[i] + 1; 1897 v = a->a + diag[i] + 1; 1898 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1899 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1900 } 1901 } 1902 if (xb == b) { 1903 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1904 } else { 1905 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1906 } 1907 } 1908 } 1909 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1910 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1911 PetscFunctionReturn(0); 1912 } 1913 1914 1915 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1916 { 1917 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1918 1919 PetscFunctionBegin; 1920 info->block_size = 1.0; 1921 info->nz_allocated = (double)a->maxnz; 1922 info->nz_used = (double)a->nz; 1923 info->nz_unneeded = (double)(a->maxnz - a->nz); 1924 info->assemblies = (double)A->num_ass; 1925 info->mallocs = (double)A->info.mallocs; 1926 info->memory = ((PetscObject)A)->mem; 1927 if (A->factortype) { 1928 info->fill_ratio_given = A->info.fill_ratio_given; 1929 info->fill_ratio_needed = A->info.fill_ratio_needed; 1930 info->factor_mallocs = A->info.factor_mallocs; 1931 } else { 1932 info->fill_ratio_given = 0; 1933 info->fill_ratio_needed = 0; 1934 info->factor_mallocs = 0; 1935 } 1936 PetscFunctionReturn(0); 1937 } 1938 1939 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1940 { 1941 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1942 PetscInt i,m = A->rmap->n - 1; 1943 PetscErrorCode ierr; 1944 const PetscScalar *xx; 1945 PetscScalar *bb; 1946 PetscInt d = 0; 1947 1948 PetscFunctionBegin; 1949 if (x && b) { 1950 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1951 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1952 for (i=0; i<N; i++) { 1953 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1954 if (rows[i] >= A->cmap->n) continue; 1955 bb[rows[i]] = diag*xx[rows[i]]; 1956 } 1957 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1958 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1959 } 1960 1961 if (a->keepnonzeropattern) { 1962 for (i=0; i<N; i++) { 1963 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1964 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1965 } 1966 if (diag != 0.0) { 1967 for (i=0; i<N; i++) { 1968 d = rows[i]; 1969 if (rows[i] >= A->cmap->n) continue; 1970 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); 1971 } 1972 for (i=0; i<N; i++) { 1973 if (rows[i] >= A->cmap->n) continue; 1974 a->a[a->diag[rows[i]]] = diag; 1975 } 1976 } 1977 } else { 1978 if (diag != 0.0) { 1979 for (i=0; i<N; i++) { 1980 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1981 if (a->ilen[rows[i]] > 0) { 1982 if (rows[i] >= A->cmap->n) { 1983 a->ilen[rows[i]] = 0; 1984 } else { 1985 a->ilen[rows[i]] = 1; 1986 a->a[a->i[rows[i]]] = diag; 1987 a->j[a->i[rows[i]]] = rows[i]; 1988 } 1989 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 1990 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1991 } 1992 } 1993 } else { 1994 for (i=0; i<N; i++) { 1995 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1996 a->ilen[rows[i]] = 0; 1997 } 1998 } 1999 A->nonzerostate++; 2000 } 2001 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2002 PetscFunctionReturn(0); 2003 } 2004 2005 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2006 { 2007 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2008 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2009 PetscErrorCode ierr; 2010 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2011 const PetscScalar *xx; 2012 PetscScalar *bb; 2013 2014 PetscFunctionBegin; 2015 if (x && b) { 2016 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2017 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2018 vecs = PETSC_TRUE; 2019 } 2020 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2021 for (i=0; i<N; i++) { 2022 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2023 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2024 2025 zeroed[rows[i]] = PETSC_TRUE; 2026 } 2027 for (i=0; i<A->rmap->n; i++) { 2028 if (!zeroed[i]) { 2029 for (j=a->i[i]; j<a->i[i+1]; j++) { 2030 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2031 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2032 a->a[j] = 0.0; 2033 } 2034 } 2035 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2036 } 2037 if (x && b) { 2038 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2039 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2040 } 2041 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2042 if (diag != 0.0) { 2043 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2044 if (missing) { 2045 for (i=0; i<N; i++) { 2046 if (rows[i] >= A->cmap->N) continue; 2047 if (a->nonew && rows[i] >= d) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D (%D)",d,rows[i]); 2048 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 2049 } 2050 } else { 2051 for (i=0; i<N; i++) { 2052 a->a[a->diag[rows[i]]] = diag; 2053 } 2054 } 2055 } 2056 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2057 PetscFunctionReturn(0); 2058 } 2059 2060 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2061 { 2062 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2063 PetscInt *itmp; 2064 2065 PetscFunctionBegin; 2066 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2067 2068 *nz = a->i[row+1] - a->i[row]; 2069 if (v) *v = a->a + a->i[row]; 2070 if (idx) { 2071 itmp = a->j + a->i[row]; 2072 if (*nz) *idx = itmp; 2073 else *idx = 0; 2074 } 2075 PetscFunctionReturn(0); 2076 } 2077 2078 /* remove this function? */ 2079 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2080 { 2081 PetscFunctionBegin; 2082 PetscFunctionReturn(0); 2083 } 2084 2085 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2086 { 2087 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2088 MatScalar *v = a->a; 2089 PetscReal sum = 0.0; 2090 PetscErrorCode ierr; 2091 PetscInt i,j; 2092 2093 PetscFunctionBegin; 2094 if (type == NORM_FROBENIUS) { 2095 #if defined(PETSC_USE_REAL___FP16) 2096 PetscBLASInt one = 1,nz = a->nz; 2097 *nrm = BLASnrm2_(&nz,v,&one); 2098 #else 2099 for (i=0; i<a->nz; i++) { 2100 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2101 } 2102 *nrm = PetscSqrtReal(sum); 2103 #endif 2104 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2105 } else if (type == NORM_1) { 2106 PetscReal *tmp; 2107 PetscInt *jj = a->j; 2108 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2109 *nrm = 0.0; 2110 for (j=0; j<a->nz; j++) { 2111 tmp[*jj++] += PetscAbsScalar(*v); v++; 2112 } 2113 for (j=0; j<A->cmap->n; j++) { 2114 if (tmp[j] > *nrm) *nrm = tmp[j]; 2115 } 2116 ierr = PetscFree(tmp);CHKERRQ(ierr); 2117 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2118 } else if (type == NORM_INFINITY) { 2119 *nrm = 0.0; 2120 for (j=0; j<A->rmap->n; j++) { 2121 v = a->a + a->i[j]; 2122 sum = 0.0; 2123 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2124 sum += PetscAbsScalar(*v); v++; 2125 } 2126 if (sum > *nrm) *nrm = sum; 2127 } 2128 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2129 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2130 PetscFunctionReturn(0); 2131 } 2132 2133 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2134 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2135 { 2136 PetscErrorCode ierr; 2137 PetscInt i,j,anzj; 2138 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2139 PetscInt an=A->cmap->N,am=A->rmap->N; 2140 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2141 2142 PetscFunctionBegin; 2143 /* Allocate space for symbolic transpose info and work array */ 2144 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 2145 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2146 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2147 2148 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2149 /* Note: offset by 1 for fast conversion into csr format. */ 2150 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2151 /* Form ati for csr format of A^T. */ 2152 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2153 2154 /* Copy ati into atfill so we have locations of the next free space in atj */ 2155 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2156 2157 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2158 for (i=0;i<am;i++) { 2159 anzj = ai[i+1] - ai[i]; 2160 for (j=0;j<anzj;j++) { 2161 atj[atfill[*aj]] = i; 2162 atfill[*aj++] += 1; 2163 } 2164 } 2165 2166 /* Clean up temporary space and complete requests. */ 2167 ierr = PetscFree(atfill);CHKERRQ(ierr); 2168 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2169 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2170 2171 b = (Mat_SeqAIJ*)((*B)->data); 2172 b->free_a = PETSC_FALSE; 2173 b->free_ij = PETSC_TRUE; 2174 b->nonew = 0; 2175 PetscFunctionReturn(0); 2176 } 2177 2178 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2179 { 2180 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2181 Mat C; 2182 PetscErrorCode ierr; 2183 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2184 MatScalar *array = a->a; 2185 2186 PetscFunctionBegin; 2187 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 2188 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2189 2190 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2191 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2192 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2193 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2194 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2195 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2196 ierr = PetscFree(col);CHKERRQ(ierr); 2197 } else { 2198 C = *B; 2199 } 2200 2201 for (i=0; i<m; i++) { 2202 len = ai[i+1]-ai[i]; 2203 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2204 array += len; 2205 aj += len; 2206 } 2207 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2208 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2209 2210 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 2211 *B = C; 2212 } else { 2213 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2214 } 2215 PetscFunctionReturn(0); 2216 } 2217 2218 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2219 { 2220 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2221 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2222 MatScalar *va,*vb; 2223 PetscErrorCode ierr; 2224 PetscInt ma,na,mb,nb, i; 2225 2226 PetscFunctionBegin; 2227 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2228 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2229 if (ma!=nb || na!=mb) { 2230 *f = PETSC_FALSE; 2231 PetscFunctionReturn(0); 2232 } 2233 aii = aij->i; bii = bij->i; 2234 adx = aij->j; bdx = bij->j; 2235 va = aij->a; vb = bij->a; 2236 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2237 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2238 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2239 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2240 2241 *f = PETSC_TRUE; 2242 for (i=0; i<ma; i++) { 2243 while (aptr[i]<aii[i+1]) { 2244 PetscInt idc,idr; 2245 PetscScalar vc,vr; 2246 /* column/row index/value */ 2247 idc = adx[aptr[i]]; 2248 idr = bdx[bptr[idc]]; 2249 vc = va[aptr[i]]; 2250 vr = vb[bptr[idc]]; 2251 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2252 *f = PETSC_FALSE; 2253 goto done; 2254 } else { 2255 aptr[i]++; 2256 if (B || i!=idc) bptr[idc]++; 2257 } 2258 } 2259 } 2260 done: 2261 ierr = PetscFree(aptr);CHKERRQ(ierr); 2262 ierr = PetscFree(bptr);CHKERRQ(ierr); 2263 PetscFunctionReturn(0); 2264 } 2265 2266 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2267 { 2268 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2269 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2270 MatScalar *va,*vb; 2271 PetscErrorCode ierr; 2272 PetscInt ma,na,mb,nb, i; 2273 2274 PetscFunctionBegin; 2275 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2276 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2277 if (ma!=nb || na!=mb) { 2278 *f = PETSC_FALSE; 2279 PetscFunctionReturn(0); 2280 } 2281 aii = aij->i; bii = bij->i; 2282 adx = aij->j; bdx = bij->j; 2283 va = aij->a; vb = bij->a; 2284 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2285 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2286 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2287 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2288 2289 *f = PETSC_TRUE; 2290 for (i=0; i<ma; i++) { 2291 while (aptr[i]<aii[i+1]) { 2292 PetscInt idc,idr; 2293 PetscScalar vc,vr; 2294 /* column/row index/value */ 2295 idc = adx[aptr[i]]; 2296 idr = bdx[bptr[idc]]; 2297 vc = va[aptr[i]]; 2298 vr = vb[bptr[idc]]; 2299 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2300 *f = PETSC_FALSE; 2301 goto done; 2302 } else { 2303 aptr[i]++; 2304 if (B || i!=idc) bptr[idc]++; 2305 } 2306 } 2307 } 2308 done: 2309 ierr = PetscFree(aptr);CHKERRQ(ierr); 2310 ierr = PetscFree(bptr);CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 2314 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2315 { 2316 PetscErrorCode ierr; 2317 2318 PetscFunctionBegin; 2319 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2324 { 2325 PetscErrorCode ierr; 2326 2327 PetscFunctionBegin; 2328 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2329 PetscFunctionReturn(0); 2330 } 2331 2332 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2333 { 2334 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2335 const PetscScalar *l,*r; 2336 PetscScalar x; 2337 MatScalar *v; 2338 PetscErrorCode ierr; 2339 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2340 const PetscInt *jj; 2341 2342 PetscFunctionBegin; 2343 if (ll) { 2344 /* The local size is used so that VecMPI can be passed to this routine 2345 by MatDiagonalScale_MPIAIJ */ 2346 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2347 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2348 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2349 v = a->a; 2350 for (i=0; i<m; i++) { 2351 x = l[i]; 2352 M = a->i[i+1] - a->i[i]; 2353 for (j=0; j<M; j++) (*v++) *= x; 2354 } 2355 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2356 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2357 } 2358 if (rr) { 2359 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2360 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2361 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2362 v = a->a; jj = a->j; 2363 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2364 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2365 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2366 } 2367 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2368 PetscFunctionReturn(0); 2369 } 2370 2371 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2372 { 2373 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2374 PetscErrorCode ierr; 2375 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2376 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2377 const PetscInt *irow,*icol; 2378 PetscInt nrows,ncols; 2379 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2380 MatScalar *a_new,*mat_a; 2381 Mat C; 2382 PetscBool stride; 2383 2384 PetscFunctionBegin; 2385 2386 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2387 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2388 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2389 2390 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2391 if (stride) { 2392 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2393 } else { 2394 first = 0; 2395 step = 0; 2396 } 2397 if (stride && step == 1) { 2398 /* special case of contiguous rows */ 2399 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2400 /* loop over new rows determining lens and starting points */ 2401 for (i=0; i<nrows; i++) { 2402 kstart = ai[irow[i]]; 2403 kend = kstart + ailen[irow[i]]; 2404 starts[i] = kstart; 2405 for (k=kstart; k<kend; k++) { 2406 if (aj[k] >= first) { 2407 starts[i] = k; 2408 break; 2409 } 2410 } 2411 sum = 0; 2412 while (k < kend) { 2413 if (aj[k++] >= first+ncols) break; 2414 sum++; 2415 } 2416 lens[i] = sum; 2417 } 2418 /* create submatrix */ 2419 if (scall == MAT_REUSE_MATRIX) { 2420 PetscInt n_cols,n_rows; 2421 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2422 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2423 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2424 C = *B; 2425 } else { 2426 PetscInt rbs,cbs; 2427 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2428 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2429 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2430 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2431 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2432 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2433 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2434 } 2435 c = (Mat_SeqAIJ*)C->data; 2436 2437 /* loop over rows inserting into submatrix */ 2438 a_new = c->a; 2439 j_new = c->j; 2440 i_new = c->i; 2441 2442 for (i=0; i<nrows; i++) { 2443 ii = starts[i]; 2444 lensi = lens[i]; 2445 for (k=0; k<lensi; k++) { 2446 *j_new++ = aj[ii+k] - first; 2447 } 2448 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2449 a_new += lensi; 2450 i_new[i+1] = i_new[i] + lensi; 2451 c->ilen[i] = lensi; 2452 } 2453 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2454 } else { 2455 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2456 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2457 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2458 for (i=0; i<ncols; i++) { 2459 #if defined(PETSC_USE_DEBUG) 2460 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); 2461 #endif 2462 smap[icol[i]] = i+1; 2463 } 2464 2465 /* determine lens of each row */ 2466 for (i=0; i<nrows; i++) { 2467 kstart = ai[irow[i]]; 2468 kend = kstart + a->ilen[irow[i]]; 2469 lens[i] = 0; 2470 for (k=kstart; k<kend; k++) { 2471 if (smap[aj[k]]) { 2472 lens[i]++; 2473 } 2474 } 2475 } 2476 /* Create and fill new matrix */ 2477 if (scall == MAT_REUSE_MATRIX) { 2478 PetscBool equal; 2479 2480 c = (Mat_SeqAIJ*)((*B)->data); 2481 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2482 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2483 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2484 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2485 C = *B; 2486 } else { 2487 PetscInt rbs,cbs; 2488 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2489 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2490 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2491 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2492 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2493 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2494 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2495 } 2496 c = (Mat_SeqAIJ*)(C->data); 2497 for (i=0; i<nrows; i++) { 2498 row = irow[i]; 2499 kstart = ai[row]; 2500 kend = kstart + a->ilen[row]; 2501 mat_i = c->i[i]; 2502 mat_j = c->j + mat_i; 2503 mat_a = c->a + mat_i; 2504 mat_ilen = c->ilen + i; 2505 for (k=kstart; k<kend; k++) { 2506 if ((tcol=smap[a->j[k]])) { 2507 *mat_j++ = tcol - 1; 2508 *mat_a++ = a->a[k]; 2509 (*mat_ilen)++; 2510 2511 } 2512 } 2513 } 2514 /* Free work space */ 2515 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2516 ierr = PetscFree(smap);CHKERRQ(ierr); 2517 ierr = PetscFree(lens);CHKERRQ(ierr); 2518 /* sort */ 2519 for (i = 0; i < nrows; i++) { 2520 PetscInt ilen; 2521 2522 mat_i = c->i[i]; 2523 mat_j = c->j + mat_i; 2524 mat_a = c->a + mat_i; 2525 ilen = c->ilen[i]; 2526 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2527 } 2528 } 2529 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2530 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2531 2532 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2533 *B = C; 2534 PetscFunctionReturn(0); 2535 } 2536 2537 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2538 { 2539 PetscErrorCode ierr; 2540 Mat B; 2541 2542 PetscFunctionBegin; 2543 if (scall == MAT_INITIAL_MATRIX) { 2544 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2545 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2546 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2547 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2548 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2549 *subMat = B; 2550 } else { 2551 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2552 } 2553 PetscFunctionReturn(0); 2554 } 2555 2556 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2557 { 2558 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2559 PetscErrorCode ierr; 2560 Mat outA; 2561 PetscBool row_identity,col_identity; 2562 2563 PetscFunctionBegin; 2564 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2565 2566 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2567 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2568 2569 outA = inA; 2570 outA->factortype = MAT_FACTOR_LU; 2571 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2572 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2573 2574 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2575 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2576 2577 a->row = row; 2578 2579 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2580 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2581 2582 a->col = col; 2583 2584 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2585 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2586 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2587 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2588 2589 if (!a->solve_work) { /* this matrix may have been factored before */ 2590 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2591 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2592 } 2593 2594 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2595 if (row_identity && col_identity) { 2596 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2597 } else { 2598 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2599 } 2600 PetscFunctionReturn(0); 2601 } 2602 2603 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2604 { 2605 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2606 PetscScalar oalpha = alpha; 2607 PetscErrorCode ierr; 2608 PetscBLASInt one = 1,bnz; 2609 2610 PetscFunctionBegin; 2611 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2612 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2613 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2614 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2615 PetscFunctionReturn(0); 2616 } 2617 2618 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2619 { 2620 PetscErrorCode ierr; 2621 PetscInt i; 2622 2623 PetscFunctionBegin; 2624 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2625 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2626 2627 for (i=0; i<submatj->nrqr; ++i) { 2628 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2629 } 2630 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2631 2632 if (submatj->rbuf1) { 2633 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2634 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2635 } 2636 2637 for (i=0; i<submatj->nrqs; ++i) { 2638 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2639 } 2640 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2641 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2642 } 2643 2644 #if defined(PETSC_USE_CTABLE) 2645 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2646 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2647 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2648 #else 2649 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2650 #endif 2651 2652 if (!submatj->allcolumns) { 2653 #if defined(PETSC_USE_CTABLE) 2654 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2655 #else 2656 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2657 #endif 2658 } 2659 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2660 2661 ierr = PetscFree(submatj);CHKERRQ(ierr); 2662 PetscFunctionReturn(0); 2663 } 2664 2665 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2666 { 2667 PetscErrorCode ierr; 2668 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2669 Mat_SubSppt *submatj = c->submatis1; 2670 2671 PetscFunctionBegin; 2672 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2673 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2674 PetscFunctionReturn(0); 2675 } 2676 2677 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2678 { 2679 PetscErrorCode ierr; 2680 PetscInt i; 2681 Mat C; 2682 Mat_SeqAIJ *c; 2683 Mat_SubSppt *submatj; 2684 2685 PetscFunctionBegin; 2686 for (i=0; i<n; i++) { 2687 C = (*mat)[i]; 2688 c = (Mat_SeqAIJ*)C->data; 2689 submatj = c->submatis1; 2690 if (submatj) { 2691 if (--((PetscObject)C)->refct <= 0) { 2692 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2693 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2694 ierr = PetscFree(C->defaultvectype);CHKERRQ(ierr); 2695 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2696 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2697 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2698 } 2699 } else { 2700 ierr = MatDestroy(&C);CHKERRQ(ierr); 2701 } 2702 } 2703 2704 /* Destroy Dummy submatrices created for reuse */ 2705 ierr = MatDestroySubMatrices_Dummy(n,mat);CHKERRQ(ierr); 2706 2707 ierr = PetscFree(*mat);CHKERRQ(ierr); 2708 PetscFunctionReturn(0); 2709 } 2710 2711 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2712 { 2713 PetscErrorCode ierr; 2714 PetscInt i; 2715 2716 PetscFunctionBegin; 2717 if (scall == MAT_INITIAL_MATRIX) { 2718 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2719 } 2720 2721 for (i=0; i<n; i++) { 2722 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2723 } 2724 PetscFunctionReturn(0); 2725 } 2726 2727 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2728 { 2729 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2730 PetscErrorCode ierr; 2731 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2732 const PetscInt *idx; 2733 PetscInt start,end,*ai,*aj; 2734 PetscBT table; 2735 2736 PetscFunctionBegin; 2737 m = A->rmap->n; 2738 ai = a->i; 2739 aj = a->j; 2740 2741 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2742 2743 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2744 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2745 2746 for (i=0; i<is_max; i++) { 2747 /* Initialize the two local arrays */ 2748 isz = 0; 2749 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2750 2751 /* Extract the indices, assume there can be duplicate entries */ 2752 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2753 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2754 2755 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2756 for (j=0; j<n; ++j) { 2757 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2758 } 2759 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2760 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2761 2762 k = 0; 2763 for (j=0; j<ov; j++) { /* for each overlap */ 2764 n = isz; 2765 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2766 row = nidx[k]; 2767 start = ai[row]; 2768 end = ai[row+1]; 2769 for (l = start; l<end; l++) { 2770 val = aj[l]; 2771 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2772 } 2773 } 2774 } 2775 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2776 } 2777 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2778 ierr = PetscFree(nidx);CHKERRQ(ierr); 2779 PetscFunctionReturn(0); 2780 } 2781 2782 /* -------------------------------------------------------------- */ 2783 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2784 { 2785 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2786 PetscErrorCode ierr; 2787 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2788 const PetscInt *row,*col; 2789 PetscInt *cnew,j,*lens; 2790 IS icolp,irowp; 2791 PetscInt *cwork = NULL; 2792 PetscScalar *vwork = NULL; 2793 2794 PetscFunctionBegin; 2795 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2796 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2797 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2798 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2799 2800 /* determine lengths of permuted rows */ 2801 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2802 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2803 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2804 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2805 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2806 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2807 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2808 ierr = PetscFree(lens);CHKERRQ(ierr); 2809 2810 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2811 for (i=0; i<m; i++) { 2812 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2813 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2814 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2815 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2816 } 2817 ierr = PetscFree(cnew);CHKERRQ(ierr); 2818 2819 (*B)->assembled = PETSC_FALSE; 2820 2821 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2822 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2823 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2824 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2825 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2826 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2827 PetscFunctionReturn(0); 2828 } 2829 2830 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2831 { 2832 PetscErrorCode ierr; 2833 2834 PetscFunctionBegin; 2835 /* If the two matrices have the same copy implementation, use fast copy. */ 2836 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2837 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2838 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2839 2840 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"); 2841 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2842 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2843 } else { 2844 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2845 } 2846 PetscFunctionReturn(0); 2847 } 2848 2849 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2850 { 2851 PetscErrorCode ierr; 2852 2853 PetscFunctionBegin; 2854 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2855 PetscFunctionReturn(0); 2856 } 2857 2858 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2859 { 2860 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2861 2862 PetscFunctionBegin; 2863 *array = a->a; 2864 PetscFunctionReturn(0); 2865 } 2866 2867 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2868 { 2869 PetscFunctionBegin; 2870 PetscFunctionReturn(0); 2871 } 2872 2873 /* 2874 Computes the number of nonzeros per row needed for preallocation when X and Y 2875 have different nonzero structure. 2876 */ 2877 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2878 { 2879 PetscInt i,j,k,nzx,nzy; 2880 2881 PetscFunctionBegin; 2882 /* Set the number of nonzeros in the new matrix */ 2883 for (i=0; i<m; i++) { 2884 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2885 nzx = xi[i+1] - xi[i]; 2886 nzy = yi[i+1] - yi[i]; 2887 nnz[i] = 0; 2888 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2889 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2890 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2891 nnz[i]++; 2892 } 2893 for (; k<nzy; k++) nnz[i]++; 2894 } 2895 PetscFunctionReturn(0); 2896 } 2897 2898 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2899 { 2900 PetscInt m = Y->rmap->N; 2901 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2902 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2903 PetscErrorCode ierr; 2904 2905 PetscFunctionBegin; 2906 /* Set the number of nonzeros in the new matrix */ 2907 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2908 PetscFunctionReturn(0); 2909 } 2910 2911 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2912 { 2913 PetscErrorCode ierr; 2914 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2915 PetscBLASInt one=1,bnz; 2916 2917 PetscFunctionBegin; 2918 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2919 if (str == SAME_NONZERO_PATTERN) { 2920 PetscScalar alpha = a; 2921 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2922 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2923 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2924 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2925 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2926 } else { 2927 Mat B; 2928 PetscInt *nnz; 2929 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2930 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2931 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2932 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2933 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2934 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2935 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2936 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2937 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2938 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2939 ierr = PetscFree(nnz);CHKERRQ(ierr); 2940 } 2941 PetscFunctionReturn(0); 2942 } 2943 2944 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2945 { 2946 #if defined(PETSC_USE_COMPLEX) 2947 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2948 PetscInt i,nz; 2949 PetscScalar *a; 2950 2951 PetscFunctionBegin; 2952 nz = aij->nz; 2953 a = aij->a; 2954 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2955 #else 2956 PetscFunctionBegin; 2957 #endif 2958 PetscFunctionReturn(0); 2959 } 2960 2961 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2962 { 2963 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2964 PetscErrorCode ierr; 2965 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2966 PetscReal atmp; 2967 PetscScalar *x; 2968 MatScalar *aa; 2969 2970 PetscFunctionBegin; 2971 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2972 aa = a->a; 2973 ai = a->i; 2974 aj = a->j; 2975 2976 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2977 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2978 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2979 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2980 for (i=0; i<m; i++) { 2981 ncols = ai[1] - ai[0]; ai++; 2982 x[i] = 0.0; 2983 for (j=0; j<ncols; j++) { 2984 atmp = PetscAbsScalar(*aa); 2985 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2986 aa++; aj++; 2987 } 2988 } 2989 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2990 PetscFunctionReturn(0); 2991 } 2992 2993 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2994 { 2995 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2996 PetscErrorCode ierr; 2997 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2998 PetscScalar *x; 2999 MatScalar *aa; 3000 3001 PetscFunctionBegin; 3002 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3003 aa = a->a; 3004 ai = a->i; 3005 aj = a->j; 3006 3007 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3008 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3009 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3010 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3011 for (i=0; i<m; i++) { 3012 ncols = ai[1] - ai[0]; ai++; 3013 if (ncols == A->cmap->n) { /* row is dense */ 3014 x[i] = *aa; if (idx) idx[i] = 0; 3015 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3016 x[i] = 0.0; 3017 if (idx) { 3018 idx[i] = 0; /* in case ncols is zero */ 3019 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 3020 if (aj[j] > j) { 3021 idx[i] = j; 3022 break; 3023 } 3024 } 3025 } 3026 } 3027 for (j=0; j<ncols; j++) { 3028 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3029 aa++; aj++; 3030 } 3031 } 3032 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3033 PetscFunctionReturn(0); 3034 } 3035 3036 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3037 { 3038 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3039 PetscErrorCode ierr; 3040 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3041 PetscReal atmp; 3042 PetscScalar *x; 3043 MatScalar *aa; 3044 3045 PetscFunctionBegin; 3046 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3047 aa = a->a; 3048 ai = a->i; 3049 aj = a->j; 3050 3051 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3052 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3053 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3054 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); 3055 for (i=0; i<m; i++) { 3056 ncols = ai[1] - ai[0]; ai++; 3057 if (ncols) { 3058 /* Get first nonzero */ 3059 for (j = 0; j < ncols; j++) { 3060 atmp = PetscAbsScalar(aa[j]); 3061 if (atmp > 1.0e-12) { 3062 x[i] = atmp; 3063 if (idx) idx[i] = aj[j]; 3064 break; 3065 } 3066 } 3067 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3068 } else { 3069 x[i] = 0.0; if (idx) idx[i] = 0; 3070 } 3071 for (j = 0; j < ncols; j++) { 3072 atmp = PetscAbsScalar(*aa); 3073 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3074 aa++; aj++; 3075 } 3076 } 3077 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3078 PetscFunctionReturn(0); 3079 } 3080 3081 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3082 { 3083 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3084 PetscErrorCode ierr; 3085 PetscInt i,j,m = A->rmap->n,ncols,n; 3086 const PetscInt *ai,*aj; 3087 PetscScalar *x; 3088 const MatScalar *aa; 3089 3090 PetscFunctionBegin; 3091 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3092 aa = a->a; 3093 ai = a->i; 3094 aj = a->j; 3095 3096 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3097 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3098 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3099 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3100 for (i=0; i<m; i++) { 3101 ncols = ai[1] - ai[0]; ai++; 3102 if (ncols == A->cmap->n) { /* row is dense */ 3103 x[i] = *aa; if (idx) idx[i] = 0; 3104 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3105 x[i] = 0.0; 3106 if (idx) { /* find first implicit 0.0 in the row */ 3107 idx[i] = 0; /* in case ncols is zero */ 3108 for (j=0; j<ncols; j++) { 3109 if (aj[j] > j) { 3110 idx[i] = j; 3111 break; 3112 } 3113 } 3114 } 3115 } 3116 for (j=0; j<ncols; j++) { 3117 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3118 aa++; aj++; 3119 } 3120 } 3121 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3122 PetscFunctionReturn(0); 3123 } 3124 3125 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3126 { 3127 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3128 PetscErrorCode ierr; 3129 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3130 MatScalar *diag,work[25],*v_work; 3131 const PetscReal shift = 0.0; 3132 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3133 3134 PetscFunctionBegin; 3135 allowzeropivot = PetscNot(A->erroriffailure); 3136 if (a->ibdiagvalid) { 3137 if (values) *values = a->ibdiag; 3138 PetscFunctionReturn(0); 3139 } 3140 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3141 if (!a->ibdiag) { 3142 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3143 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3144 } 3145 diag = a->ibdiag; 3146 if (values) *values = a->ibdiag; 3147 /* factor and invert each block */ 3148 switch (bs) { 3149 case 1: 3150 for (i=0; i<mbs; i++) { 3151 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3152 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3153 if (allowzeropivot) { 3154 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3155 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3156 A->factorerror_zeropivot_row = i; 3157 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3158 } 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); 3159 } 3160 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3161 } 3162 break; 3163 case 2: 3164 for (i=0; i<mbs; i++) { 3165 ij[0] = 2*i; ij[1] = 2*i + 1; 3166 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3167 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3168 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3169 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3170 diag += 4; 3171 } 3172 break; 3173 case 3: 3174 for (i=0; i<mbs; i++) { 3175 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3176 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3177 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3178 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3179 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3180 diag += 9; 3181 } 3182 break; 3183 case 4: 3184 for (i=0; i<mbs; i++) { 3185 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3186 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3187 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3188 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3189 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3190 diag += 16; 3191 } 3192 break; 3193 case 5: 3194 for (i=0; i<mbs; i++) { 3195 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3196 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3197 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3198 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3199 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3200 diag += 25; 3201 } 3202 break; 3203 case 6: 3204 for (i=0; i<mbs; i++) { 3205 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; 3206 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3207 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3208 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3209 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3210 diag += 36; 3211 } 3212 break; 3213 case 7: 3214 for (i=0; i<mbs; i++) { 3215 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; 3216 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3217 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3218 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3219 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3220 diag += 49; 3221 } 3222 break; 3223 default: 3224 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3225 for (i=0; i<mbs; i++) { 3226 for (j=0; j<bs; j++) { 3227 IJ[j] = bs*i + j; 3228 } 3229 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3230 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3231 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3232 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3233 diag += bs2; 3234 } 3235 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3236 } 3237 a->ibdiagvalid = PETSC_TRUE; 3238 PetscFunctionReturn(0); 3239 } 3240 3241 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3242 { 3243 PetscErrorCode ierr; 3244 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3245 PetscScalar a; 3246 PetscInt m,n,i,j,col; 3247 3248 PetscFunctionBegin; 3249 if (!x->assembled) { 3250 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3251 for (i=0; i<m; i++) { 3252 for (j=0; j<aij->imax[i]; j++) { 3253 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3254 col = (PetscInt)(n*PetscRealPart(a)); 3255 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3256 } 3257 } 3258 } else { 3259 for (i=0; i<aij->nz; i++) {ierr = PetscRandomGetValue(rctx,aij->a+i);CHKERRQ(ierr);} 3260 } 3261 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3262 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3263 PetscFunctionReturn(0); 3264 } 3265 3266 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3267 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3268 { 3269 PetscErrorCode ierr; 3270 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3271 PetscScalar a; 3272 PetscInt m,n,i,j,col,nskip; 3273 3274 PetscFunctionBegin; 3275 nskip = high - low; 3276 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3277 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3278 for (i=0; i<m; i++) { 3279 for (j=0; j<aij->imax[i]; j++) { 3280 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3281 col = (PetscInt)(n*PetscRealPart(a)); 3282 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3283 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3284 } 3285 } 3286 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3287 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3288 PetscFunctionReturn(0); 3289 } 3290 3291 3292 /* -------------------------------------------------------------------*/ 3293 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3294 MatGetRow_SeqAIJ, 3295 MatRestoreRow_SeqAIJ, 3296 MatMult_SeqAIJ, 3297 /* 4*/ MatMultAdd_SeqAIJ, 3298 MatMultTranspose_SeqAIJ, 3299 MatMultTransposeAdd_SeqAIJ, 3300 0, 3301 0, 3302 0, 3303 /* 10*/ 0, 3304 MatLUFactor_SeqAIJ, 3305 0, 3306 MatSOR_SeqAIJ, 3307 MatTranspose_SeqAIJ_FAST, 3308 /*1 5*/ MatGetInfo_SeqAIJ, 3309 MatEqual_SeqAIJ, 3310 MatGetDiagonal_SeqAIJ, 3311 MatDiagonalScale_SeqAIJ, 3312 MatNorm_SeqAIJ, 3313 /* 20*/ 0, 3314 MatAssemblyEnd_SeqAIJ, 3315 MatSetOption_SeqAIJ, 3316 MatZeroEntries_SeqAIJ, 3317 /* 24*/ MatZeroRows_SeqAIJ, 3318 0, 3319 0, 3320 0, 3321 0, 3322 /* 29*/ MatSetUp_SeqAIJ, 3323 0, 3324 0, 3325 0, 3326 0, 3327 /* 34*/ MatDuplicate_SeqAIJ, 3328 0, 3329 0, 3330 MatILUFactor_SeqAIJ, 3331 0, 3332 /* 39*/ MatAXPY_SeqAIJ, 3333 MatCreateSubMatrices_SeqAIJ, 3334 MatIncreaseOverlap_SeqAIJ, 3335 MatGetValues_SeqAIJ, 3336 MatCopy_SeqAIJ, 3337 /* 44*/ MatGetRowMax_SeqAIJ, 3338 MatScale_SeqAIJ, 3339 MatShift_SeqAIJ, 3340 MatDiagonalSet_SeqAIJ, 3341 MatZeroRowsColumns_SeqAIJ, 3342 /* 49*/ MatSetRandom_SeqAIJ, 3343 MatGetRowIJ_SeqAIJ, 3344 MatRestoreRowIJ_SeqAIJ, 3345 MatGetColumnIJ_SeqAIJ, 3346 MatRestoreColumnIJ_SeqAIJ, 3347 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3348 0, 3349 0, 3350 MatPermute_SeqAIJ, 3351 0, 3352 /* 59*/ 0, 3353 MatDestroy_SeqAIJ, 3354 MatView_SeqAIJ, 3355 0, 3356 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3357 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3358 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3359 0, 3360 0, 3361 0, 3362 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3363 MatGetRowMinAbs_SeqAIJ, 3364 0, 3365 0, 3366 0, 3367 /* 74*/ 0, 3368 MatFDColoringApply_AIJ, 3369 0, 3370 0, 3371 0, 3372 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3373 0, 3374 0, 3375 0, 3376 MatLoad_SeqAIJ, 3377 /* 84*/ MatIsSymmetric_SeqAIJ, 3378 MatIsHermitian_SeqAIJ, 3379 0, 3380 0, 3381 0, 3382 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3383 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3384 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3385 MatPtAP_SeqAIJ_SeqAIJ, 3386 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3387 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3388 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3389 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3390 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3391 0, 3392 /* 99*/ 0, 3393 0, 3394 0, 3395 MatConjugate_SeqAIJ, 3396 0, 3397 /*104*/ MatSetValuesRow_SeqAIJ, 3398 MatRealPart_SeqAIJ, 3399 MatImaginaryPart_SeqAIJ, 3400 0, 3401 0, 3402 /*109*/ MatMatSolve_SeqAIJ, 3403 0, 3404 MatGetRowMin_SeqAIJ, 3405 0, 3406 MatMissingDiagonal_SeqAIJ, 3407 /*114*/ 0, 3408 0, 3409 0, 3410 0, 3411 0, 3412 /*119*/ 0, 3413 0, 3414 0, 3415 0, 3416 MatGetMultiProcBlock_SeqAIJ, 3417 /*124*/ MatFindNonzeroRows_SeqAIJ, 3418 MatGetColumnNorms_SeqAIJ, 3419 MatInvertBlockDiagonal_SeqAIJ, 3420 MatInvertVariableBlockDiagonal_SeqAIJ, 3421 0, 3422 /*129*/ 0, 3423 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3424 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3425 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3426 MatTransposeColoringCreate_SeqAIJ, 3427 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3428 MatTransColoringApplyDenToSp_SeqAIJ, 3429 MatRARt_SeqAIJ_SeqAIJ, 3430 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3431 MatRARtNumeric_SeqAIJ_SeqAIJ, 3432 /*139*/0, 3433 0, 3434 0, 3435 MatFDColoringSetUp_SeqXAIJ, 3436 MatFindOffBlockDiagonalEntries_SeqAIJ, 3437 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3438 MatDestroySubMatrices_SeqAIJ 3439 }; 3440 3441 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3442 { 3443 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3444 PetscInt i,nz,n; 3445 3446 PetscFunctionBegin; 3447 nz = aij->maxnz; 3448 n = mat->rmap->n; 3449 for (i=0; i<nz; i++) { 3450 aij->j[i] = indices[i]; 3451 } 3452 aij->nz = nz; 3453 for (i=0; i<n; i++) { 3454 aij->ilen[i] = aij->imax[i]; 3455 } 3456 PetscFunctionReturn(0); 3457 } 3458 3459 /* 3460 * When a sparse matrix has many zero columns, we should compact them out to save the space 3461 * This happens in MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3462 * */ 3463 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3464 { 3465 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3466 PetscTable gid1_lid1; 3467 PetscTablePosition tpos; 3468 PetscInt gid,lid,i,j,ncols,ec; 3469 PetscInt *garray; 3470 PetscErrorCode ierr; 3471 3472 PetscFunctionBegin; 3473 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3474 PetscValidPointer(mapping,2); 3475 /* use a table */ 3476 ierr = PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr); 3477 ec = 0; 3478 for (i=0; i<mat->rmap->n; i++) { 3479 ncols = aij->i[i+1] - aij->i[i]; 3480 for (j=0; j<ncols; j++) { 3481 PetscInt data,gid1 = aij->j[aij->i[i] + j] + 1; 3482 ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr); 3483 if (!data) { 3484 /* one based table */ 3485 ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr); 3486 } 3487 } 3488 } 3489 /* form array of columns we need */ 3490 ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr); 3491 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 3492 while (tpos) { 3493 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr); 3494 gid--; 3495 lid--; 3496 garray[lid] = gid; 3497 } 3498 ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */ 3499 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 3500 for (i=0; i<ec; i++) { 3501 ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 3502 } 3503 /* compact out the extra columns in B */ 3504 for (i=0; i<mat->rmap->n; i++) { 3505 ncols = aij->i[i+1] - aij->i[i]; 3506 for (j=0; j<ncols; j++) { 3507 PetscInt gid1 = aij->j[aij->i[i] + j] + 1; 3508 ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr); 3509 lid--; 3510 aij->j[aij->i[i] + j] = lid; 3511 } 3512 } 3513 mat->cmap->n = mat->cmap->N = ec; 3514 mat->cmap->bs = 1; 3515 3516 ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr); 3517 ierr = PetscLayoutSetUp((mat->cmap));CHKERRQ(ierr); 3518 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);CHKERRQ(ierr); 3519 ierr = ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);CHKERRQ(ierr); 3520 PetscFunctionReturn(0); 3521 } 3522 3523 /*@ 3524 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3525 in the matrix. 3526 3527 Input Parameters: 3528 + mat - the SeqAIJ matrix 3529 - indices - the column indices 3530 3531 Level: advanced 3532 3533 Notes: 3534 This can be called if you have precomputed the nonzero structure of the 3535 matrix and want to provide it to the matrix object to improve the performance 3536 of the MatSetValues() operation. 3537 3538 You MUST have set the correct numbers of nonzeros per row in the call to 3539 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3540 3541 MUST be called before any calls to MatSetValues(); 3542 3543 The indices should start with zero, not one. 3544 3545 @*/ 3546 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3547 { 3548 PetscErrorCode ierr; 3549 3550 PetscFunctionBegin; 3551 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3552 PetscValidPointer(indices,2); 3553 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3554 PetscFunctionReturn(0); 3555 } 3556 3557 /* ----------------------------------------------------------------------------------------*/ 3558 3559 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3560 { 3561 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3562 PetscErrorCode ierr; 3563 size_t nz = aij->i[mat->rmap->n]; 3564 3565 PetscFunctionBegin; 3566 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3567 3568 /* allocate space for values if not already there */ 3569 if (!aij->saved_values) { 3570 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3571 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3572 } 3573 3574 /* copy values over */ 3575 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3576 PetscFunctionReturn(0); 3577 } 3578 3579 /*@ 3580 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3581 example, reuse of the linear part of a Jacobian, while recomputing the 3582 nonlinear portion. 3583 3584 Collect on Mat 3585 3586 Input Parameters: 3587 . mat - the matrix (currently only AIJ matrices support this option) 3588 3589 Level: advanced 3590 3591 Common Usage, with SNESSolve(): 3592 $ Create Jacobian matrix 3593 $ Set linear terms into matrix 3594 $ Apply boundary conditions to matrix, at this time matrix must have 3595 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3596 $ boundary conditions again will not change the nonzero structure 3597 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3598 $ ierr = MatStoreValues(mat); 3599 $ Call SNESSetJacobian() with matrix 3600 $ In your Jacobian routine 3601 $ ierr = MatRetrieveValues(mat); 3602 $ Set nonlinear terms in matrix 3603 3604 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3605 $ // build linear portion of Jacobian 3606 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3607 $ ierr = MatStoreValues(mat); 3608 $ loop over nonlinear iterations 3609 $ ierr = MatRetrieveValues(mat); 3610 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3611 $ // call MatAssemblyBegin/End() on matrix 3612 $ Solve linear system with Jacobian 3613 $ endloop 3614 3615 Notes: 3616 Matrix must already be assemblied before calling this routine 3617 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3618 calling this routine. 3619 3620 When this is called multiple times it overwrites the previous set of stored values 3621 and does not allocated additional space. 3622 3623 .seealso: MatRetrieveValues() 3624 3625 @*/ 3626 PetscErrorCode MatStoreValues(Mat mat) 3627 { 3628 PetscErrorCode ierr; 3629 3630 PetscFunctionBegin; 3631 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3632 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3633 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3634 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3635 PetscFunctionReturn(0); 3636 } 3637 3638 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3639 { 3640 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3641 PetscErrorCode ierr; 3642 PetscInt nz = aij->i[mat->rmap->n]; 3643 3644 PetscFunctionBegin; 3645 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3646 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3647 /* copy values over */ 3648 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3649 PetscFunctionReturn(0); 3650 } 3651 3652 /*@ 3653 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3654 example, reuse of the linear part of a Jacobian, while recomputing the 3655 nonlinear portion. 3656 3657 Collect on Mat 3658 3659 Input Parameters: 3660 . mat - the matrix (currently only AIJ matrices support this option) 3661 3662 Level: advanced 3663 3664 .seealso: MatStoreValues() 3665 3666 @*/ 3667 PetscErrorCode MatRetrieveValues(Mat mat) 3668 { 3669 PetscErrorCode ierr; 3670 3671 PetscFunctionBegin; 3672 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3673 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3674 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3675 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3676 PetscFunctionReturn(0); 3677 } 3678 3679 3680 /* --------------------------------------------------------------------------------*/ 3681 /*@C 3682 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3683 (the default parallel PETSc format). For good matrix assembly performance 3684 the user should preallocate the matrix storage by setting the parameter nz 3685 (or the array nnz). By setting these parameters accurately, performance 3686 during matrix assembly can be increased by more than a factor of 50. 3687 3688 Collective on MPI_Comm 3689 3690 Input Parameters: 3691 + comm - MPI communicator, set to PETSC_COMM_SELF 3692 . m - number of rows 3693 . n - number of columns 3694 . nz - number of nonzeros per row (same for all rows) 3695 - nnz - array containing the number of nonzeros in the various rows 3696 (possibly different for each row) or NULL 3697 3698 Output Parameter: 3699 . A - the matrix 3700 3701 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3702 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3703 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3704 3705 Notes: 3706 If nnz is given then nz is ignored 3707 3708 The AIJ format (also called the Yale sparse matrix format or 3709 compressed row storage), is fully compatible with standard Fortran 77 3710 storage. That is, the stored row and column indices can begin at 3711 either one (as in Fortran) or zero. See the users' manual for details. 3712 3713 Specify the preallocated storage with either nz or nnz (not both). 3714 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3715 allocation. For large problems you MUST preallocate memory or you 3716 will get TERRIBLE performance, see the users' manual chapter on matrices. 3717 3718 By default, this format uses inodes (identical nodes) when possible, to 3719 improve numerical efficiency of matrix-vector products and solves. We 3720 search for consecutive rows with the same nonzero structure, thereby 3721 reusing matrix information to achieve increased efficiency. 3722 3723 Options Database Keys: 3724 + -mat_no_inode - Do not use inodes 3725 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3726 3727 Level: intermediate 3728 3729 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3730 3731 @*/ 3732 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3733 { 3734 PetscErrorCode ierr; 3735 3736 PetscFunctionBegin; 3737 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3738 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3739 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3740 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3741 PetscFunctionReturn(0); 3742 } 3743 3744 /*@C 3745 MatSeqAIJSetPreallocation - For good matrix assembly performance 3746 the user should preallocate the matrix storage by setting the parameter nz 3747 (or the array nnz). By setting these parameters accurately, performance 3748 during matrix assembly can be increased by more than a factor of 50. 3749 3750 Collective on MPI_Comm 3751 3752 Input Parameters: 3753 + B - The matrix 3754 . nz - number of nonzeros per row (same for all rows) 3755 - nnz - array containing the number of nonzeros in the various rows 3756 (possibly different for each row) or NULL 3757 3758 Notes: 3759 If nnz is given then nz is ignored 3760 3761 The AIJ format (also called the Yale sparse matrix format or 3762 compressed row storage), is fully compatible with standard Fortran 77 3763 storage. That is, the stored row and column indices can begin at 3764 either one (as in Fortran) or zero. See the users' manual for details. 3765 3766 Specify the preallocated storage with either nz or nnz (not both). 3767 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3768 allocation. For large problems you MUST preallocate memory or you 3769 will get TERRIBLE performance, see the users' manual chapter on matrices. 3770 3771 You can call MatGetInfo() to get information on how effective the preallocation was; 3772 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3773 You can also run with the option -info and look for messages with the string 3774 malloc in them to see if additional memory allocation was needed. 3775 3776 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3777 entries or columns indices 3778 3779 By default, this format uses inodes (identical nodes) when possible, to 3780 improve numerical efficiency of matrix-vector products and solves. We 3781 search for consecutive rows with the same nonzero structure, thereby 3782 reusing matrix information to achieve increased efficiency. 3783 3784 Options Database Keys: 3785 + -mat_no_inode - Do not use inodes 3786 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3787 3788 Level: intermediate 3789 3790 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3791 3792 @*/ 3793 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3794 { 3795 PetscErrorCode ierr; 3796 3797 PetscFunctionBegin; 3798 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3799 PetscValidType(B,1); 3800 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3801 PetscFunctionReturn(0); 3802 } 3803 3804 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3805 { 3806 Mat_SeqAIJ *b; 3807 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3808 PetscErrorCode ierr; 3809 PetscInt i; 3810 3811 PetscFunctionBegin; 3812 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3813 if (nz == MAT_SKIP_ALLOCATION) { 3814 skipallocation = PETSC_TRUE; 3815 nz = 0; 3816 } 3817 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3818 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3819 3820 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3821 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3822 if (nnz) { 3823 for (i=0; i<B->rmap->n; i++) { 3824 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]); 3825 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); 3826 } 3827 } 3828 3829 B->preallocated = PETSC_TRUE; 3830 3831 b = (Mat_SeqAIJ*)B->data; 3832 3833 if (!skipallocation) { 3834 if (!b->imax) { 3835 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3836 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3837 } 3838 if (!b->ipre) { 3839 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3840 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3841 } 3842 if (!nnz) { 3843 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3844 else if (nz < 0) nz = 1; 3845 nz = PetscMin(nz,B->cmap->n); 3846 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3847 nz = nz*B->rmap->n; 3848 } else { 3849 nz = 0; 3850 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3851 } 3852 /* b->ilen will count nonzeros in each row so far. */ 3853 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3854 3855 /* allocate the matrix space */ 3856 /* FIXME: should B's old memory be unlogged? */ 3857 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3858 if (B->structure_only) { 3859 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3860 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3861 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3862 } else { 3863 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3864 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3865 } 3866 b->i[0] = 0; 3867 for (i=1; i<B->rmap->n+1; i++) { 3868 b->i[i] = b->i[i-1] + b->imax[i-1]; 3869 } 3870 if (B->structure_only) { 3871 b->singlemalloc = PETSC_FALSE; 3872 b->free_a = PETSC_FALSE; 3873 } else { 3874 b->singlemalloc = PETSC_TRUE; 3875 b->free_a = PETSC_TRUE; 3876 } 3877 b->free_ij = PETSC_TRUE; 3878 } else { 3879 b->free_a = PETSC_FALSE; 3880 b->free_ij = PETSC_FALSE; 3881 } 3882 3883 if (b->ipre && nnz != b->ipre && b->imax) { 3884 /* reserve user-requested sparsity */ 3885 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3886 } 3887 3888 3889 b->nz = 0; 3890 b->maxnz = nz; 3891 B->info.nz_unneeded = (double)b->maxnz; 3892 if (realalloc) { 3893 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3894 } 3895 B->was_assembled = PETSC_FALSE; 3896 B->assembled = PETSC_FALSE; 3897 PetscFunctionReturn(0); 3898 } 3899 3900 3901 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3902 { 3903 Mat_SeqAIJ *a; 3904 PetscInt i; 3905 PetscErrorCode ierr; 3906 3907 PetscFunctionBegin; 3908 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3909 a = (Mat_SeqAIJ*)A->data; 3910 /* if no saved info, we error out */ 3911 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3912 3913 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"); 3914 3915 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3916 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3917 a->i[0] = 0; 3918 for (i=1; i<A->rmap->n+1; i++) { 3919 a->i[i] = a->i[i-1] + a->imax[i-1]; 3920 } 3921 A->preallocated = PETSC_TRUE; 3922 a->nz = 0; 3923 a->maxnz = a->i[A->rmap->n]; 3924 A->info.nz_unneeded = (double)a->maxnz; 3925 A->was_assembled = PETSC_FALSE; 3926 A->assembled = PETSC_FALSE; 3927 PetscFunctionReturn(0); 3928 } 3929 3930 /*@ 3931 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3932 3933 Input Parameters: 3934 + B - the matrix 3935 . i - the indices into j for the start of each row (starts with zero) 3936 . j - the column indices for each row (starts with zero) these must be sorted for each row 3937 - v - optional values in the matrix 3938 3939 Level: developer 3940 3941 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3942 3943 .keywords: matrix, aij, compressed row, sparse, sequential 3944 3945 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ 3946 @*/ 3947 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3948 { 3949 PetscErrorCode ierr; 3950 3951 PetscFunctionBegin; 3952 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3953 PetscValidType(B,1); 3954 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3955 PetscFunctionReturn(0); 3956 } 3957 3958 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3959 { 3960 PetscInt i; 3961 PetscInt m,n; 3962 PetscInt nz; 3963 PetscInt *nnz, nz_max = 0; 3964 PetscScalar *values; 3965 PetscErrorCode ierr; 3966 3967 PetscFunctionBegin; 3968 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3969 3970 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3971 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3972 3973 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3974 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3975 for (i = 0; i < m; i++) { 3976 nz = Ii[i+1]- Ii[i]; 3977 nz_max = PetscMax(nz_max, nz); 3978 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3979 nnz[i] = nz; 3980 } 3981 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3982 ierr = PetscFree(nnz);CHKERRQ(ierr); 3983 3984 if (v) { 3985 values = (PetscScalar*) v; 3986 } else { 3987 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3988 } 3989 3990 for (i = 0; i < m; i++) { 3991 nz = Ii[i+1] - Ii[i]; 3992 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3993 } 3994 3995 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3996 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3997 3998 if (!v) { 3999 ierr = PetscFree(values);CHKERRQ(ierr); 4000 } 4001 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 4002 PetscFunctionReturn(0); 4003 } 4004 4005 #include <../src/mat/impls/dense/seq/dense.h> 4006 #include <petsc/private/kernels/petscaxpy.h> 4007 4008 /* 4009 Computes (B'*A')' since computing B*A directly is untenable 4010 4011 n p p 4012 ( ) ( ) ( ) 4013 m ( A ) * n ( B ) = m ( C ) 4014 ( ) ( ) ( ) 4015 4016 */ 4017 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4018 { 4019 PetscErrorCode ierr; 4020 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4021 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4022 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4023 PetscInt i,n,m,q,p; 4024 const PetscInt *ii,*idx; 4025 const PetscScalar *b,*a,*a_q; 4026 PetscScalar *c,*c_q; 4027 4028 PetscFunctionBegin; 4029 m = A->rmap->n; 4030 n = A->cmap->n; 4031 p = B->cmap->n; 4032 a = sub_a->v; 4033 b = sub_b->a; 4034 c = sub_c->v; 4035 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 4036 4037 ii = sub_b->i; 4038 idx = sub_b->j; 4039 for (i=0; i<n; i++) { 4040 q = ii[i+1] - ii[i]; 4041 while (q-->0) { 4042 c_q = c + m*(*idx); 4043 a_q = a + m*i; 4044 PetscKernelAXPY(c_q,*b,a_q,m); 4045 idx++; 4046 b++; 4047 } 4048 } 4049 PetscFunctionReturn(0); 4050 } 4051 4052 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 4053 { 4054 PetscErrorCode ierr; 4055 PetscInt m=A->rmap->n,n=B->cmap->n; 4056 Mat Cmat; 4057 4058 PetscFunctionBegin; 4059 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); 4060 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 4061 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 4062 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 4063 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 4064 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 4065 4066 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4067 4068 *C = Cmat; 4069 PetscFunctionReturn(0); 4070 } 4071 4072 /* ----------------------------------------------------------------*/ 4073 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 4074 { 4075 PetscErrorCode ierr; 4076 4077 PetscFunctionBegin; 4078 if (scall == MAT_INITIAL_MATRIX) { 4079 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4080 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 4081 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4082 } 4083 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4084 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 4085 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4086 PetscFunctionReturn(0); 4087 } 4088 4089 4090 /*MC 4091 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4092 based on compressed sparse row format. 4093 4094 Options Database Keys: 4095 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4096 4097 Level: beginner 4098 4099 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 4100 M*/ 4101 4102 /*MC 4103 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4104 4105 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4106 and MATMPIAIJ otherwise. As a result, for single process communicators, 4107 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 4108 for communicators controlling multiple processes. It is recommended that you call both of 4109 the above preallocation routines for simplicity. 4110 4111 Options Database Keys: 4112 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4113 4114 Developer Notes: 4115 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4116 enough exist. 4117 4118 Level: beginner 4119 4120 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 4121 M*/ 4122 4123 /*MC 4124 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4125 4126 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4127 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4128 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4129 for communicators controlling multiple processes. It is recommended that you call both of 4130 the above preallocation routines for simplicity. 4131 4132 Options Database Keys: 4133 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4134 4135 Level: beginner 4136 4137 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4138 M*/ 4139 4140 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4141 #if defined(PETSC_HAVE_ELEMENTAL) 4142 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4143 #endif 4144 #if defined(PETSC_HAVE_HYPRE) 4145 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4146 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 4147 #endif 4148 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 4149 4150 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4151 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4152 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 4153 4154 /*@C 4155 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 4156 4157 Not Collective 4158 4159 Input Parameter: 4160 . mat - a MATSEQAIJ matrix 4161 4162 Output Parameter: 4163 . array - pointer to the data 4164 4165 Level: intermediate 4166 4167 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4168 @*/ 4169 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4170 { 4171 PetscErrorCode ierr; 4172 4173 PetscFunctionBegin; 4174 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4175 PetscFunctionReturn(0); 4176 } 4177 4178 /*@C 4179 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4180 4181 Not Collective 4182 4183 Input Parameter: 4184 . mat - a MATSEQAIJ matrix 4185 4186 Output Parameter: 4187 . nz - the maximum number of nonzeros in any row 4188 4189 Level: intermediate 4190 4191 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4192 @*/ 4193 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4194 { 4195 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4196 4197 PetscFunctionBegin; 4198 *nz = aij->rmax; 4199 PetscFunctionReturn(0); 4200 } 4201 4202 /*@C 4203 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4204 4205 Not Collective 4206 4207 Input Parameters: 4208 . mat - a MATSEQAIJ matrix 4209 . array - pointer to the data 4210 4211 Level: intermediate 4212 4213 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4214 @*/ 4215 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4216 { 4217 PetscErrorCode ierr; 4218 4219 PetscFunctionBegin; 4220 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4221 PetscFunctionReturn(0); 4222 } 4223 4224 #if defined(PETSC_HAVE_CUDA) 4225 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat); 4226 #endif 4227 4228 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4229 { 4230 Mat_SeqAIJ *b; 4231 PetscErrorCode ierr; 4232 PetscMPIInt size; 4233 4234 PetscFunctionBegin; 4235 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4236 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4237 4238 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4239 4240 B->data = (void*)b; 4241 4242 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4243 4244 b->row = 0; 4245 b->col = 0; 4246 b->icol = 0; 4247 b->reallocs = 0; 4248 b->ignorezeroentries = PETSC_FALSE; 4249 b->roworiented = PETSC_TRUE; 4250 b->nonew = 0; 4251 b->diag = 0; 4252 b->solve_work = 0; 4253 B->spptr = 0; 4254 b->saved_values = 0; 4255 b->idiag = 0; 4256 b->mdiag = 0; 4257 b->ssor_work = 0; 4258 b->omega = 1.0; 4259 b->fshift = 0.0; 4260 b->idiagvalid = PETSC_FALSE; 4261 b->ibdiagvalid = PETSC_FALSE; 4262 b->keepnonzeropattern = PETSC_FALSE; 4263 4264 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4265 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4266 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4267 4268 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4269 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4270 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4271 #endif 4272 4273 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4274 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4275 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4276 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4277 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4278 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4279 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4280 #if defined(PETSC_HAVE_MKL_SPARSE) 4281 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4282 #endif 4283 #if defined(PETSC_HAVE_CUDA) 4284 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);CHKERRQ(ierr); 4285 #endif 4286 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4287 #if defined(PETSC_HAVE_ELEMENTAL) 4288 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4289 #endif 4290 #if defined(PETSC_HAVE_HYPRE) 4291 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4292 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4293 #endif 4294 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4295 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4296 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4297 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4298 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4299 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4300 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4301 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4302 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4303 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4304 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4305 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4306 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaij_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 4307 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4308 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4309 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4310 PetscFunctionReturn(0); 4311 } 4312 4313 /* 4314 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4315 */ 4316 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4317 { 4318 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4319 PetscErrorCode ierr; 4320 PetscInt i,m = A->rmap->n; 4321 4322 PetscFunctionBegin; 4323 c = (Mat_SeqAIJ*)C->data; 4324 4325 C->factortype = A->factortype; 4326 c->row = 0; 4327 c->col = 0; 4328 c->icol = 0; 4329 c->reallocs = 0; 4330 4331 C->assembled = PETSC_TRUE; 4332 4333 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4334 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4335 4336 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4337 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4338 for (i=0; i<m; i++) { 4339 c->imax[i] = a->imax[i]; 4340 c->ilen[i] = a->ilen[i]; 4341 } 4342 4343 /* allocate the matrix space */ 4344 if (mallocmatspace) { 4345 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4346 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4347 4348 c->singlemalloc = PETSC_TRUE; 4349 4350 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4351 if (m > 0) { 4352 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4353 if (cpvalues == MAT_COPY_VALUES) { 4354 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4355 } else { 4356 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4357 } 4358 } 4359 } 4360 4361 c->ignorezeroentries = a->ignorezeroentries; 4362 c->roworiented = a->roworiented; 4363 c->nonew = a->nonew; 4364 if (a->diag) { 4365 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4366 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4367 for (i=0; i<m; i++) { 4368 c->diag[i] = a->diag[i]; 4369 } 4370 } else c->diag = 0; 4371 4372 c->solve_work = 0; 4373 c->saved_values = 0; 4374 c->idiag = 0; 4375 c->ssor_work = 0; 4376 c->keepnonzeropattern = a->keepnonzeropattern; 4377 c->free_a = PETSC_TRUE; 4378 c->free_ij = PETSC_TRUE; 4379 4380 c->rmax = a->rmax; 4381 c->nz = a->nz; 4382 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4383 C->preallocated = PETSC_TRUE; 4384 4385 c->compressedrow.use = a->compressedrow.use; 4386 c->compressedrow.nrows = a->compressedrow.nrows; 4387 if (a->compressedrow.use) { 4388 i = a->compressedrow.nrows; 4389 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4390 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4391 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4392 } else { 4393 c->compressedrow.use = PETSC_FALSE; 4394 c->compressedrow.i = NULL; 4395 c->compressedrow.rindex = NULL; 4396 } 4397 c->nonzerorowcnt = a->nonzerorowcnt; 4398 C->nonzerostate = A->nonzerostate; 4399 4400 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4401 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4402 PetscFunctionReturn(0); 4403 } 4404 4405 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4406 { 4407 PetscErrorCode ierr; 4408 4409 PetscFunctionBegin; 4410 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4411 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4412 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4413 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4414 } 4415 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4416 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4417 PetscFunctionReturn(0); 4418 } 4419 4420 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4421 { 4422 PetscBool isbinary, ishdf5; 4423 PetscErrorCode ierr; 4424 4425 PetscFunctionBegin; 4426 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4427 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4428 /* force binary viewer to load .info file if it has not yet done so */ 4429 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4430 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 4431 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 4432 if (isbinary) { 4433 ierr = MatLoad_SeqAIJ_Binary(newMat,viewer);CHKERRQ(ierr); 4434 } else if (ishdf5) { 4435 #if defined(PETSC_HAVE_HDF5) 4436 ierr = MatLoad_AIJ_HDF5(newMat,viewer);CHKERRQ(ierr); 4437 #else 4438 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4439 #endif 4440 } else { 4441 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); 4442 } 4443 PetscFunctionReturn(0); 4444 } 4445 4446 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat newMat, PetscViewer viewer) 4447 { 4448 Mat_SeqAIJ *a; 4449 PetscErrorCode ierr; 4450 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4451 int fd; 4452 PetscMPIInt size; 4453 MPI_Comm comm; 4454 PetscInt bs = newMat->rmap->bs; 4455 4456 PetscFunctionBegin; 4457 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4458 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4459 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4460 4461 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4462 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4463 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4464 if (bs < 0) bs = 1; 4465 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4466 4467 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4468 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4469 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4470 M = header[1]; N = header[2]; nz = header[3]; 4471 4472 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4473 4474 /* read in row lengths */ 4475 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4476 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4477 4478 /* check if sum of rowlengths is same as nz */ 4479 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4480 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); 4481 4482 /* set global size if not set already*/ 4483 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4484 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4485 } else { 4486 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4487 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4488 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4489 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4490 } 4491 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); 4492 } 4493 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4494 a = (Mat_SeqAIJ*)newMat->data; 4495 4496 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4497 4498 /* read in nonzero values */ 4499 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4500 4501 /* set matrix "i" values */ 4502 a->i[0] = 0; 4503 for (i=1; i<= M; i++) { 4504 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4505 a->ilen[i-1] = rowlengths[i-1]; 4506 } 4507 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4508 4509 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4510 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4511 PetscFunctionReturn(0); 4512 } 4513 4514 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4515 { 4516 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4517 PetscErrorCode ierr; 4518 #if defined(PETSC_USE_COMPLEX) 4519 PetscInt k; 4520 #endif 4521 4522 PetscFunctionBegin; 4523 /* If the matrix dimensions are not equal,or no of nonzeros */ 4524 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4525 *flg = PETSC_FALSE; 4526 PetscFunctionReturn(0); 4527 } 4528 4529 /* if the a->i are the same */ 4530 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4531 if (!*flg) PetscFunctionReturn(0); 4532 4533 /* if a->j are the same */ 4534 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4535 if (!*flg) PetscFunctionReturn(0); 4536 4537 /* if a->a are the same */ 4538 #if defined(PETSC_USE_COMPLEX) 4539 for (k=0; k<a->nz; k++) { 4540 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4541 *flg = PETSC_FALSE; 4542 PetscFunctionReturn(0); 4543 } 4544 } 4545 #else 4546 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4547 #endif 4548 PetscFunctionReturn(0); 4549 } 4550 4551 /*@ 4552 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4553 provided by the user. 4554 4555 Collective on MPI_Comm 4556 4557 Input Parameters: 4558 + comm - must be an MPI communicator of size 1 4559 . m - number of rows 4560 . n - number of columns 4561 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 4562 . j - column indices 4563 - a - matrix values 4564 4565 Output Parameter: 4566 . mat - the matrix 4567 4568 Level: intermediate 4569 4570 Notes: 4571 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4572 once the matrix is destroyed and not before 4573 4574 You cannot set new nonzero locations into this matrix, that will generate an error. 4575 4576 The i and j indices are 0 based 4577 4578 The format which is used for the sparse matrix input, is equivalent to a 4579 row-major ordering.. i.e for the following matrix, the input data expected is 4580 as shown 4581 4582 $ 1 0 0 4583 $ 2 0 3 4584 $ 4 5 6 4585 $ 4586 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4587 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4588 $ v = {1,2,3,4,5,6} [size = 6] 4589 4590 4591 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4592 4593 @*/ 4594 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4595 { 4596 PetscErrorCode ierr; 4597 PetscInt ii; 4598 Mat_SeqAIJ *aij; 4599 #if defined(PETSC_USE_DEBUG) 4600 PetscInt jj; 4601 #endif 4602 4603 PetscFunctionBegin; 4604 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4605 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4606 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4607 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4608 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4609 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4610 aij = (Mat_SeqAIJ*)(*mat)->data; 4611 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4612 4613 aij->i = i; 4614 aij->j = j; 4615 aij->a = a; 4616 aij->singlemalloc = PETSC_FALSE; 4617 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4618 aij->free_a = PETSC_FALSE; 4619 aij->free_ij = PETSC_FALSE; 4620 4621 for (ii=0; ii<m; ii++) { 4622 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4623 #if defined(PETSC_USE_DEBUG) 4624 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]); 4625 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4626 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); 4627 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); 4628 } 4629 #endif 4630 } 4631 #if defined(PETSC_USE_DEBUG) 4632 for (ii=0; ii<aij->i[m]; ii++) { 4633 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4634 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]); 4635 } 4636 #endif 4637 4638 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4639 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4640 PetscFunctionReturn(0); 4641 } 4642 /*@C 4643 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4644 provided by the user. 4645 4646 Collective on MPI_Comm 4647 4648 Input Parameters: 4649 + comm - must be an MPI communicator of size 1 4650 . m - number of rows 4651 . n - number of columns 4652 . i - row indices 4653 . j - column indices 4654 . a - matrix values 4655 . nz - number of nonzeros 4656 - idx - 0 or 1 based 4657 4658 Output Parameter: 4659 . mat - the matrix 4660 4661 Level: intermediate 4662 4663 Notes: 4664 The i and j indices are 0 based 4665 4666 The format which is used for the sparse matrix input, is equivalent to a 4667 row-major ordering.. i.e for the following matrix, the input data expected is 4668 as shown: 4669 4670 1 0 0 4671 2 0 3 4672 4 5 6 4673 4674 i = {0,1,1,2,2,2} 4675 j = {0,0,2,0,1,2} 4676 v = {1,2,3,4,5,6} 4677 4678 4679 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4680 4681 @*/ 4682 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4683 { 4684 PetscErrorCode ierr; 4685 PetscInt ii, *nnz, one = 1,row,col; 4686 4687 4688 PetscFunctionBegin; 4689 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4690 for (ii = 0; ii < nz; ii++) { 4691 nnz[i[ii] - !!idx] += 1; 4692 } 4693 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4694 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4695 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4696 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4697 for (ii = 0; ii < nz; ii++) { 4698 if (idx) { 4699 row = i[ii] - 1; 4700 col = j[ii] - 1; 4701 } else { 4702 row = i[ii]; 4703 col = j[ii]; 4704 } 4705 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4706 } 4707 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4708 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4709 ierr = PetscFree(nnz);CHKERRQ(ierr); 4710 PetscFunctionReturn(0); 4711 } 4712 4713 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4714 { 4715 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4716 PetscErrorCode ierr; 4717 4718 PetscFunctionBegin; 4719 a->idiagvalid = PETSC_FALSE; 4720 a->ibdiagvalid = PETSC_FALSE; 4721 4722 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4723 PetscFunctionReturn(0); 4724 } 4725 4726 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4727 { 4728 PetscErrorCode ierr; 4729 PetscMPIInt size; 4730 4731 PetscFunctionBegin; 4732 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4733 if (size == 1) { 4734 if (scall == MAT_INITIAL_MATRIX) { 4735 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4736 } else { 4737 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4738 } 4739 } else { 4740 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4741 } 4742 PetscFunctionReturn(0); 4743 } 4744 4745 /* 4746 Permute A into C's *local* index space using rowemb,colemb. 4747 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4748 of [0,m), colemb is in [0,n). 4749 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4750 */ 4751 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4752 { 4753 /* If making this function public, change the error returned in this function away from _PLIB. */ 4754 PetscErrorCode ierr; 4755 Mat_SeqAIJ *Baij; 4756 PetscBool seqaij; 4757 PetscInt m,n,*nz,i,j,count; 4758 PetscScalar v; 4759 const PetscInt *rowindices,*colindices; 4760 4761 PetscFunctionBegin; 4762 if (!B) PetscFunctionReturn(0); 4763 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4764 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4765 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4766 if (rowemb) { 4767 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4768 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); 4769 } else { 4770 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4771 } 4772 if (colemb) { 4773 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4774 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); 4775 } else { 4776 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4777 } 4778 4779 Baij = (Mat_SeqAIJ*)(B->data); 4780 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4781 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4782 for (i=0; i<B->rmap->n; i++) { 4783 nz[i] = Baij->i[i+1] - Baij->i[i]; 4784 } 4785 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4786 ierr = PetscFree(nz);CHKERRQ(ierr); 4787 } 4788 if (pattern == SUBSET_NONZERO_PATTERN) { 4789 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4790 } 4791 count = 0; 4792 rowindices = NULL; 4793 colindices = NULL; 4794 if (rowemb) { 4795 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4796 } 4797 if (colemb) { 4798 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4799 } 4800 for (i=0; i<B->rmap->n; i++) { 4801 PetscInt row; 4802 row = i; 4803 if (rowindices) row = rowindices[i]; 4804 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4805 PetscInt col; 4806 col = Baij->j[count]; 4807 if (colindices) col = colindices[col]; 4808 v = Baij->a[count]; 4809 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4810 ++count; 4811 } 4812 } 4813 /* FIXME: set C's nonzerostate correctly. */ 4814 /* Assembly for C is necessary. */ 4815 C->preallocated = PETSC_TRUE; 4816 C->assembled = PETSC_TRUE; 4817 C->was_assembled = PETSC_FALSE; 4818 PetscFunctionReturn(0); 4819 } 4820 4821 PetscFunctionList MatSeqAIJList = NULL; 4822 4823 /*@C 4824 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4825 4826 Collective on Mat 4827 4828 Input Parameters: 4829 + mat - the matrix object 4830 - matype - matrix type 4831 4832 Options Database Key: 4833 . -mat_seqai_type <method> - for example seqaijcrl 4834 4835 4836 Level: intermediate 4837 4838 .keywords: Mat, MatType, set, method 4839 4840 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4841 @*/ 4842 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4843 { 4844 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4845 PetscBool sametype; 4846 4847 PetscFunctionBegin; 4848 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4849 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4850 if (sametype) PetscFunctionReturn(0); 4851 4852 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4853 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4854 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4855 PetscFunctionReturn(0); 4856 } 4857 4858 4859 /*@C 4860 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4861 4862 Not Collective 4863 4864 Input Parameters: 4865 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4866 - function - routine to convert to subtype 4867 4868 Notes: 4869 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4870 4871 4872 Then, your matrix can be chosen with the procedural interface at runtime via the option 4873 $ -mat_seqaij_type my_mat 4874 4875 Level: advanced 4876 4877 .keywords: Mat, register 4878 4879 .seealso: MatSeqAIJRegisterAll() 4880 4881 4882 Level: advanced 4883 @*/ 4884 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4885 { 4886 PetscErrorCode ierr; 4887 4888 PetscFunctionBegin; 4889 ierr = MatInitializePackage();CHKERRQ(ierr); 4890 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4891 PetscFunctionReturn(0); 4892 } 4893 4894 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4895 4896 /*@C 4897 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4898 4899 Not Collective 4900 4901 Level: advanced 4902 4903 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4904 4905 .keywords: KSP, register, all 4906 4907 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4908 @*/ 4909 PetscErrorCode MatSeqAIJRegisterAll(void) 4910 { 4911 PetscErrorCode ierr; 4912 4913 PetscFunctionBegin; 4914 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4915 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4916 4917 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4918 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4919 ierr = MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4920 #if defined(PETSC_HAVE_MKL_SPARSE) 4921 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4922 #endif 4923 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4924 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4925 #endif 4926 PetscFunctionReturn(0); 4927 } 4928 4929 /* 4930 Special version for direct calls from Fortran 4931 */ 4932 #include <petsc/private/fortranimpl.h> 4933 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4934 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4935 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4936 #define matsetvaluesseqaij_ matsetvaluesseqaij 4937 #endif 4938 4939 /* Change these macros so can be used in void function */ 4940 #undef CHKERRQ 4941 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4942 #undef SETERRQ2 4943 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4944 #undef SETERRQ3 4945 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4946 4947 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) 4948 { 4949 Mat A = *AA; 4950 PetscInt m = *mm, n = *nn; 4951 InsertMode is = *isis; 4952 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4953 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4954 PetscInt *imax,*ai,*ailen; 4955 PetscErrorCode ierr; 4956 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4957 MatScalar *ap,value,*aa; 4958 PetscBool ignorezeroentries = a->ignorezeroentries; 4959 PetscBool roworiented = a->roworiented; 4960 4961 PetscFunctionBegin; 4962 MatCheckPreallocated(A,1); 4963 imax = a->imax; 4964 ai = a->i; 4965 ailen = a->ilen; 4966 aj = a->j; 4967 aa = a->a; 4968 4969 for (k=0; k<m; k++) { /* loop over added rows */ 4970 row = im[k]; 4971 if (row < 0) continue; 4972 #if defined(PETSC_USE_DEBUG) 4973 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4974 #endif 4975 rp = aj + ai[row]; ap = aa + ai[row]; 4976 rmax = imax[row]; nrow = ailen[row]; 4977 low = 0; 4978 high = nrow; 4979 for (l=0; l<n; l++) { /* loop over added columns */ 4980 if (in[l] < 0) continue; 4981 #if defined(PETSC_USE_DEBUG) 4982 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4983 #endif 4984 col = in[l]; 4985 if (roworiented) value = v[l + k*n]; 4986 else value = v[k + l*m]; 4987 4988 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4989 4990 if (col <= lastcol) low = 0; 4991 else high = nrow; 4992 lastcol = col; 4993 while (high-low > 5) { 4994 t = (low+high)/2; 4995 if (rp[t] > col) high = t; 4996 else low = t; 4997 } 4998 for (i=low; i<high; i++) { 4999 if (rp[i] > col) break; 5000 if (rp[i] == col) { 5001 if (is == ADD_VALUES) ap[i] += value; 5002 else ap[i] = value; 5003 goto noinsert; 5004 } 5005 } 5006 if (value == 0.0 && ignorezeroentries) goto noinsert; 5007 if (nonew == 1) goto noinsert; 5008 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5009 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5010 N = nrow++ - 1; a->nz++; high++; 5011 /* shift up all the later entries in this row */ 5012 for (ii=N; ii>=i; ii--) { 5013 rp[ii+1] = rp[ii]; 5014 ap[ii+1] = ap[ii]; 5015 } 5016 rp[i] = col; 5017 ap[i] = value; 5018 A->nonzerostate++; 5019 noinsert:; 5020 low = i + 1; 5021 } 5022 ailen[row] = nrow; 5023 } 5024 PetscFunctionReturnVoid(); 5025 } 5026