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