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