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 = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2623 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2624 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2625 } 2626 } else { 2627 ierr = MatDestroy(&C);CHKERRQ(ierr); 2628 } 2629 } 2630 2631 ierr = PetscFree(*mat);CHKERRQ(ierr); 2632 PetscFunctionReturn(0); 2633 } 2634 2635 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2636 { 2637 PetscErrorCode ierr; 2638 PetscInt i; 2639 2640 PetscFunctionBegin; 2641 if (scall == MAT_INITIAL_MATRIX) { 2642 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2643 } 2644 2645 for (i=0; i<n; i++) { 2646 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2647 } 2648 PetscFunctionReturn(0); 2649 } 2650 2651 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2652 { 2653 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2654 PetscErrorCode ierr; 2655 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2656 const PetscInt *idx; 2657 PetscInt start,end,*ai,*aj; 2658 PetscBT table; 2659 2660 PetscFunctionBegin; 2661 m = A->rmap->n; 2662 ai = a->i; 2663 aj = a->j; 2664 2665 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2666 2667 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2668 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2669 2670 for (i=0; i<is_max; i++) { 2671 /* Initialize the two local arrays */ 2672 isz = 0; 2673 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2674 2675 /* Extract the indices, assume there can be duplicate entries */ 2676 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2677 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2678 2679 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2680 for (j=0; j<n; ++j) { 2681 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2682 } 2683 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2684 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2685 2686 k = 0; 2687 for (j=0; j<ov; j++) { /* for each overlap */ 2688 n = isz; 2689 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2690 row = nidx[k]; 2691 start = ai[row]; 2692 end = ai[row+1]; 2693 for (l = start; l<end; l++) { 2694 val = aj[l]; 2695 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2696 } 2697 } 2698 } 2699 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2700 } 2701 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2702 ierr = PetscFree(nidx);CHKERRQ(ierr); 2703 PetscFunctionReturn(0); 2704 } 2705 2706 /* -------------------------------------------------------------- */ 2707 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2708 { 2709 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2710 PetscErrorCode ierr; 2711 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2712 const PetscInt *row,*col; 2713 PetscInt *cnew,j,*lens; 2714 IS icolp,irowp; 2715 PetscInt *cwork = NULL; 2716 PetscScalar *vwork = NULL; 2717 2718 PetscFunctionBegin; 2719 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2720 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2721 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2722 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2723 2724 /* determine lengths of permuted rows */ 2725 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2726 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2727 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2728 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2729 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2730 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2731 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2732 ierr = PetscFree(lens);CHKERRQ(ierr); 2733 2734 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2735 for (i=0; i<m; i++) { 2736 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2737 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2738 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2739 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2740 } 2741 ierr = PetscFree(cnew);CHKERRQ(ierr); 2742 2743 (*B)->assembled = PETSC_FALSE; 2744 2745 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2746 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2747 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2748 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2749 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2750 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2751 PetscFunctionReturn(0); 2752 } 2753 2754 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2755 { 2756 PetscErrorCode ierr; 2757 2758 PetscFunctionBegin; 2759 /* If the two matrices have the same copy implementation, use fast copy. */ 2760 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2761 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2762 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2763 2764 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"); 2765 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2766 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2767 } else { 2768 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2769 } 2770 PetscFunctionReturn(0); 2771 } 2772 2773 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2774 { 2775 PetscErrorCode ierr; 2776 2777 PetscFunctionBegin; 2778 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2779 PetscFunctionReturn(0); 2780 } 2781 2782 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2783 { 2784 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2785 2786 PetscFunctionBegin; 2787 *array = a->a; 2788 PetscFunctionReturn(0); 2789 } 2790 2791 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2792 { 2793 PetscFunctionBegin; 2794 PetscFunctionReturn(0); 2795 } 2796 2797 /* 2798 Computes the number of nonzeros per row needed for preallocation when X and Y 2799 have different nonzero structure. 2800 */ 2801 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2802 { 2803 PetscInt i,j,k,nzx,nzy; 2804 2805 PetscFunctionBegin; 2806 /* Set the number of nonzeros in the new matrix */ 2807 for (i=0; i<m; i++) { 2808 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2809 nzx = xi[i+1] - xi[i]; 2810 nzy = yi[i+1] - yi[i]; 2811 nnz[i] = 0; 2812 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2813 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2814 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2815 nnz[i]++; 2816 } 2817 for (; k<nzy; k++) nnz[i]++; 2818 } 2819 PetscFunctionReturn(0); 2820 } 2821 2822 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2823 { 2824 PetscInt m = Y->rmap->N; 2825 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2826 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2827 PetscErrorCode ierr; 2828 2829 PetscFunctionBegin; 2830 /* Set the number of nonzeros in the new matrix */ 2831 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2832 PetscFunctionReturn(0); 2833 } 2834 2835 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2836 { 2837 PetscErrorCode ierr; 2838 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2839 PetscBLASInt one=1,bnz; 2840 2841 PetscFunctionBegin; 2842 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2843 if (str == SAME_NONZERO_PATTERN) { 2844 PetscScalar alpha = a; 2845 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2846 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2847 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2848 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2849 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2850 } else { 2851 Mat B; 2852 PetscInt *nnz; 2853 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2854 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2855 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2856 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2857 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2858 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2859 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2860 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2861 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2862 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2863 ierr = PetscFree(nnz);CHKERRQ(ierr); 2864 } 2865 PetscFunctionReturn(0); 2866 } 2867 2868 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2869 { 2870 #if defined(PETSC_USE_COMPLEX) 2871 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2872 PetscInt i,nz; 2873 PetscScalar *a; 2874 2875 PetscFunctionBegin; 2876 nz = aij->nz; 2877 a = aij->a; 2878 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2879 #else 2880 PetscFunctionBegin; 2881 #endif 2882 PetscFunctionReturn(0); 2883 } 2884 2885 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2886 { 2887 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2888 PetscErrorCode ierr; 2889 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2890 PetscReal atmp; 2891 PetscScalar *x; 2892 MatScalar *aa; 2893 2894 PetscFunctionBegin; 2895 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2896 aa = a->a; 2897 ai = a->i; 2898 aj = a->j; 2899 2900 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2901 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2902 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2903 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2904 for (i=0; i<m; i++) { 2905 ncols = ai[1] - ai[0]; ai++; 2906 x[i] = 0.0; 2907 for (j=0; j<ncols; j++) { 2908 atmp = PetscAbsScalar(*aa); 2909 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2910 aa++; aj++; 2911 } 2912 } 2913 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2914 PetscFunctionReturn(0); 2915 } 2916 2917 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2918 { 2919 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2920 PetscErrorCode ierr; 2921 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2922 PetscScalar *x; 2923 MatScalar *aa; 2924 2925 PetscFunctionBegin; 2926 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2927 aa = a->a; 2928 ai = a->i; 2929 aj = a->j; 2930 2931 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2932 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2933 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2934 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2935 for (i=0; i<m; i++) { 2936 ncols = ai[1] - ai[0]; ai++; 2937 if (ncols == A->cmap->n) { /* row is dense */ 2938 x[i] = *aa; if (idx) idx[i] = 0; 2939 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2940 x[i] = 0.0; 2941 if (idx) { 2942 idx[i] = 0; /* in case ncols is zero */ 2943 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2944 if (aj[j] > j) { 2945 idx[i] = j; 2946 break; 2947 } 2948 } 2949 } 2950 } 2951 for (j=0; j<ncols; j++) { 2952 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2953 aa++; aj++; 2954 } 2955 } 2956 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2957 PetscFunctionReturn(0); 2958 } 2959 2960 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2961 { 2962 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2963 PetscErrorCode ierr; 2964 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2965 PetscReal atmp; 2966 PetscScalar *x; 2967 MatScalar *aa; 2968 2969 PetscFunctionBegin; 2970 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2971 aa = a->a; 2972 ai = a->i; 2973 aj = a->j; 2974 2975 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2976 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2977 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2978 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); 2979 for (i=0; i<m; i++) { 2980 ncols = ai[1] - ai[0]; ai++; 2981 if (ncols) { 2982 /* Get first nonzero */ 2983 for (j = 0; j < ncols; j++) { 2984 atmp = PetscAbsScalar(aa[j]); 2985 if (atmp > 1.0e-12) { 2986 x[i] = atmp; 2987 if (idx) idx[i] = aj[j]; 2988 break; 2989 } 2990 } 2991 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2992 } else { 2993 x[i] = 0.0; if (idx) idx[i] = 0; 2994 } 2995 for (j = 0; j < ncols; j++) { 2996 atmp = PetscAbsScalar(*aa); 2997 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2998 aa++; aj++; 2999 } 3000 } 3001 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3002 PetscFunctionReturn(0); 3003 } 3004 3005 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3006 { 3007 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3008 PetscErrorCode ierr; 3009 PetscInt i,j,m = A->rmap->n,ncols,n; 3010 const PetscInt *ai,*aj; 3011 PetscScalar *x; 3012 const MatScalar *aa; 3013 3014 PetscFunctionBegin; 3015 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3016 aa = a->a; 3017 ai = a->i; 3018 aj = a->j; 3019 3020 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3021 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3022 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3023 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3024 for (i=0; i<m; i++) { 3025 ncols = ai[1] - ai[0]; ai++; 3026 if (ncols == A->cmap->n) { /* row is dense */ 3027 x[i] = *aa; if (idx) idx[i] = 0; 3028 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3029 x[i] = 0.0; 3030 if (idx) { /* find first implicit 0.0 in the row */ 3031 idx[i] = 0; /* in case ncols is zero */ 3032 for (j=0; j<ncols; j++) { 3033 if (aj[j] > j) { 3034 idx[i] = j; 3035 break; 3036 } 3037 } 3038 } 3039 } 3040 for (j=0; j<ncols; j++) { 3041 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3042 aa++; aj++; 3043 } 3044 } 3045 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3046 PetscFunctionReturn(0); 3047 } 3048 3049 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3050 { 3051 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3052 PetscErrorCode ierr; 3053 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3054 MatScalar *diag,work[25],*v_work; 3055 const PetscReal shift = 0.0; 3056 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3057 3058 PetscFunctionBegin; 3059 allowzeropivot = PetscNot(A->erroriffailure); 3060 if (a->ibdiagvalid) { 3061 if (values) *values = a->ibdiag; 3062 PetscFunctionReturn(0); 3063 } 3064 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3065 if (!a->ibdiag) { 3066 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3067 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3068 } 3069 diag = a->ibdiag; 3070 if (values) *values = a->ibdiag; 3071 /* factor and invert each block */ 3072 switch (bs) { 3073 case 1: 3074 for (i=0; i<mbs; i++) { 3075 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3076 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3077 if (allowzeropivot) { 3078 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3079 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3080 A->factorerror_zeropivot_row = i; 3081 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3082 } 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); 3083 } 3084 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3085 } 3086 break; 3087 case 2: 3088 for (i=0; i<mbs; i++) { 3089 ij[0] = 2*i; ij[1] = 2*i + 1; 3090 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3091 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3092 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3093 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3094 diag += 4; 3095 } 3096 break; 3097 case 3: 3098 for (i=0; i<mbs; i++) { 3099 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3100 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3101 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3102 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3103 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3104 diag += 9; 3105 } 3106 break; 3107 case 4: 3108 for (i=0; i<mbs; i++) { 3109 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3110 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3111 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3112 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3113 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3114 diag += 16; 3115 } 3116 break; 3117 case 5: 3118 for (i=0; i<mbs; i++) { 3119 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3120 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3121 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3122 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3123 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3124 diag += 25; 3125 } 3126 break; 3127 case 6: 3128 for (i=0; i<mbs; i++) { 3129 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; 3130 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3131 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3132 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3133 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3134 diag += 36; 3135 } 3136 break; 3137 case 7: 3138 for (i=0; i<mbs; i++) { 3139 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; 3140 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3141 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3142 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3143 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3144 diag += 49; 3145 } 3146 break; 3147 default: 3148 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3149 for (i=0; i<mbs; i++) { 3150 for (j=0; j<bs; j++) { 3151 IJ[j] = bs*i + j; 3152 } 3153 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3154 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3155 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3156 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3157 diag += bs2; 3158 } 3159 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3160 } 3161 a->ibdiagvalid = PETSC_TRUE; 3162 PetscFunctionReturn(0); 3163 } 3164 3165 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3166 { 3167 PetscErrorCode ierr; 3168 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3169 PetscScalar a; 3170 PetscInt m,n,i,j,col; 3171 3172 PetscFunctionBegin; 3173 if (!x->assembled) { 3174 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3175 for (i=0; i<m; i++) { 3176 for (j=0; j<aij->imax[i]; j++) { 3177 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3178 col = (PetscInt)(n*PetscRealPart(a)); 3179 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3180 } 3181 } 3182 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3183 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3184 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3185 PetscFunctionReturn(0); 3186 } 3187 3188 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3189 { 3190 PetscErrorCode ierr; 3191 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3192 3193 PetscFunctionBegin; 3194 if (!Y->preallocated || !aij->nz) { 3195 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3196 } 3197 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3198 PetscFunctionReturn(0); 3199 } 3200 3201 /* -------------------------------------------------------------------*/ 3202 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3203 MatGetRow_SeqAIJ, 3204 MatRestoreRow_SeqAIJ, 3205 MatMult_SeqAIJ, 3206 /* 4*/ MatMultAdd_SeqAIJ, 3207 MatMultTranspose_SeqAIJ, 3208 MatMultTransposeAdd_SeqAIJ, 3209 0, 3210 0, 3211 0, 3212 /* 10*/ 0, 3213 MatLUFactor_SeqAIJ, 3214 0, 3215 MatSOR_SeqAIJ, 3216 MatTranspose_SeqAIJ_FAST, 3217 /*1 5*/ MatGetInfo_SeqAIJ, 3218 MatEqual_SeqAIJ, 3219 MatGetDiagonal_SeqAIJ, 3220 MatDiagonalScale_SeqAIJ, 3221 MatNorm_SeqAIJ, 3222 /* 20*/ 0, 3223 MatAssemblyEnd_SeqAIJ, 3224 MatSetOption_SeqAIJ, 3225 MatZeroEntries_SeqAIJ, 3226 /* 24*/ MatZeroRows_SeqAIJ, 3227 0, 3228 0, 3229 0, 3230 0, 3231 /* 29*/ MatSetUp_SeqAIJ, 3232 0, 3233 0, 3234 0, 3235 0, 3236 /* 34*/ MatDuplicate_SeqAIJ, 3237 0, 3238 0, 3239 MatILUFactor_SeqAIJ, 3240 0, 3241 /* 39*/ MatAXPY_SeqAIJ, 3242 MatCreateSubMatrices_SeqAIJ, 3243 MatIncreaseOverlap_SeqAIJ, 3244 MatGetValues_SeqAIJ, 3245 MatCopy_SeqAIJ, 3246 /* 44*/ MatGetRowMax_SeqAIJ, 3247 MatScale_SeqAIJ, 3248 MatShift_SeqAIJ, 3249 MatDiagonalSet_SeqAIJ, 3250 MatZeroRowsColumns_SeqAIJ, 3251 /* 49*/ MatSetRandom_SeqAIJ, 3252 MatGetRowIJ_SeqAIJ, 3253 MatRestoreRowIJ_SeqAIJ, 3254 MatGetColumnIJ_SeqAIJ, 3255 MatRestoreColumnIJ_SeqAIJ, 3256 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3257 0, 3258 0, 3259 MatPermute_SeqAIJ, 3260 0, 3261 /* 59*/ 0, 3262 MatDestroy_SeqAIJ, 3263 MatView_SeqAIJ, 3264 0, 3265 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3266 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3267 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3268 0, 3269 0, 3270 0, 3271 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3272 MatGetRowMinAbs_SeqAIJ, 3273 0, 3274 0, 3275 0, 3276 /* 74*/ 0, 3277 MatFDColoringApply_AIJ, 3278 0, 3279 0, 3280 0, 3281 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3282 0, 3283 0, 3284 0, 3285 MatLoad_SeqAIJ, 3286 /* 84*/ MatIsSymmetric_SeqAIJ, 3287 MatIsHermitian_SeqAIJ, 3288 0, 3289 0, 3290 0, 3291 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3292 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3293 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3294 MatPtAP_SeqAIJ_SeqAIJ, 3295 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3296 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3297 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3298 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3299 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3300 0, 3301 /* 99*/ 0, 3302 0, 3303 0, 3304 MatConjugate_SeqAIJ, 3305 0, 3306 /*104*/ MatSetValuesRow_SeqAIJ, 3307 MatRealPart_SeqAIJ, 3308 MatImaginaryPart_SeqAIJ, 3309 0, 3310 0, 3311 /*109*/ MatMatSolve_SeqAIJ, 3312 0, 3313 MatGetRowMin_SeqAIJ, 3314 0, 3315 MatMissingDiagonal_SeqAIJ, 3316 /*114*/ 0, 3317 0, 3318 0, 3319 0, 3320 0, 3321 /*119*/ 0, 3322 0, 3323 0, 3324 0, 3325 MatGetMultiProcBlock_SeqAIJ, 3326 /*124*/ MatFindNonzeroRows_SeqAIJ, 3327 MatGetColumnNorms_SeqAIJ, 3328 MatInvertBlockDiagonal_SeqAIJ, 3329 MatInvertVariableBlockDiagonal_SeqAIJ, 3330 0, 3331 /*129*/ 0, 3332 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3333 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3334 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3335 MatTransposeColoringCreate_SeqAIJ, 3336 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3337 MatTransColoringApplyDenToSp_SeqAIJ, 3338 MatRARt_SeqAIJ_SeqAIJ, 3339 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3340 MatRARtNumeric_SeqAIJ_SeqAIJ, 3341 /*139*/0, 3342 0, 3343 0, 3344 MatFDColoringSetUp_SeqXAIJ, 3345 MatFindOffBlockDiagonalEntries_SeqAIJ, 3346 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3347 MatDestroySubMatrices_SeqAIJ 3348 }; 3349 3350 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3351 { 3352 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3353 PetscInt i,nz,n; 3354 3355 PetscFunctionBegin; 3356 nz = aij->maxnz; 3357 n = mat->rmap->n; 3358 for (i=0; i<nz; i++) { 3359 aij->j[i] = indices[i]; 3360 } 3361 aij->nz = nz; 3362 for (i=0; i<n; i++) { 3363 aij->ilen[i] = aij->imax[i]; 3364 } 3365 PetscFunctionReturn(0); 3366 } 3367 3368 /*@ 3369 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3370 in the matrix. 3371 3372 Input Parameters: 3373 + mat - the SeqAIJ matrix 3374 - indices - the column indices 3375 3376 Level: advanced 3377 3378 Notes: 3379 This can be called if you have precomputed the nonzero structure of the 3380 matrix and want to provide it to the matrix object to improve the performance 3381 of the MatSetValues() operation. 3382 3383 You MUST have set the correct numbers of nonzeros per row in the call to 3384 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3385 3386 MUST be called before any calls to MatSetValues(); 3387 3388 The indices should start with zero, not one. 3389 3390 @*/ 3391 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3392 { 3393 PetscErrorCode ierr; 3394 3395 PetscFunctionBegin; 3396 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3397 PetscValidPointer(indices,2); 3398 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3399 PetscFunctionReturn(0); 3400 } 3401 3402 /* ----------------------------------------------------------------------------------------*/ 3403 3404 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3405 { 3406 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3407 PetscErrorCode ierr; 3408 size_t nz = aij->i[mat->rmap->n]; 3409 3410 PetscFunctionBegin; 3411 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3412 3413 /* allocate space for values if not already there */ 3414 if (!aij->saved_values) { 3415 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3416 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3417 } 3418 3419 /* copy values over */ 3420 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3421 PetscFunctionReturn(0); 3422 } 3423 3424 /*@ 3425 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3426 example, reuse of the linear part of a Jacobian, while recomputing the 3427 nonlinear portion. 3428 3429 Collect on Mat 3430 3431 Input Parameters: 3432 . mat - the matrix (currently only AIJ matrices support this option) 3433 3434 Level: advanced 3435 3436 Common Usage, with SNESSolve(): 3437 $ Create Jacobian matrix 3438 $ Set linear terms into matrix 3439 $ Apply boundary conditions to matrix, at this time matrix must have 3440 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3441 $ boundary conditions again will not change the nonzero structure 3442 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3443 $ ierr = MatStoreValues(mat); 3444 $ Call SNESSetJacobian() with matrix 3445 $ In your Jacobian routine 3446 $ ierr = MatRetrieveValues(mat); 3447 $ Set nonlinear terms in matrix 3448 3449 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3450 $ // build linear portion of Jacobian 3451 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3452 $ ierr = MatStoreValues(mat); 3453 $ loop over nonlinear iterations 3454 $ ierr = MatRetrieveValues(mat); 3455 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3456 $ // call MatAssemblyBegin/End() on matrix 3457 $ Solve linear system with Jacobian 3458 $ endloop 3459 3460 Notes: 3461 Matrix must already be assemblied before calling this routine 3462 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3463 calling this routine. 3464 3465 When this is called multiple times it overwrites the previous set of stored values 3466 and does not allocated additional space. 3467 3468 .seealso: MatRetrieveValues() 3469 3470 @*/ 3471 PetscErrorCode MatStoreValues(Mat mat) 3472 { 3473 PetscErrorCode ierr; 3474 3475 PetscFunctionBegin; 3476 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3477 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3478 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3479 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3480 PetscFunctionReturn(0); 3481 } 3482 3483 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3484 { 3485 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3486 PetscErrorCode ierr; 3487 PetscInt nz = aij->i[mat->rmap->n]; 3488 3489 PetscFunctionBegin; 3490 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3491 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3492 /* copy values over */ 3493 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3494 PetscFunctionReturn(0); 3495 } 3496 3497 /*@ 3498 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3499 example, reuse of the linear part of a Jacobian, while recomputing the 3500 nonlinear portion. 3501 3502 Collect on Mat 3503 3504 Input Parameters: 3505 . mat - the matrix (currently only AIJ matrices support this option) 3506 3507 Level: advanced 3508 3509 .seealso: MatStoreValues() 3510 3511 @*/ 3512 PetscErrorCode MatRetrieveValues(Mat mat) 3513 { 3514 PetscErrorCode ierr; 3515 3516 PetscFunctionBegin; 3517 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3518 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3519 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3520 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3521 PetscFunctionReturn(0); 3522 } 3523 3524 3525 /* --------------------------------------------------------------------------------*/ 3526 /*@C 3527 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3528 (the default parallel PETSc format). For good matrix assembly performance 3529 the user should preallocate the matrix storage by setting the parameter nz 3530 (or the array nnz). By setting these parameters accurately, performance 3531 during matrix assembly can be increased by more than a factor of 50. 3532 3533 Collective on MPI_Comm 3534 3535 Input Parameters: 3536 + comm - MPI communicator, set to PETSC_COMM_SELF 3537 . m - number of rows 3538 . n - number of columns 3539 . nz - number of nonzeros per row (same for all rows) 3540 - nnz - array containing the number of nonzeros in the various rows 3541 (possibly different for each row) or NULL 3542 3543 Output Parameter: 3544 . A - the matrix 3545 3546 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3547 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3548 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3549 3550 Notes: 3551 If nnz is given then nz is ignored 3552 3553 The AIJ format (also called the Yale sparse matrix format or 3554 compressed row storage), is fully compatible with standard Fortran 77 3555 storage. That is, the stored row and column indices can begin at 3556 either one (as in Fortran) or zero. See the users' manual for details. 3557 3558 Specify the preallocated storage with either nz or nnz (not both). 3559 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3560 allocation. For large problems you MUST preallocate memory or you 3561 will get TERRIBLE performance, see the users' manual chapter on matrices. 3562 3563 By default, this format uses inodes (identical nodes) when possible, to 3564 improve numerical efficiency of matrix-vector products and solves. We 3565 search for consecutive rows with the same nonzero structure, thereby 3566 reusing matrix information to achieve increased efficiency. 3567 3568 Options Database Keys: 3569 + -mat_no_inode - Do not use inodes 3570 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3571 3572 Level: intermediate 3573 3574 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3575 3576 @*/ 3577 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3578 { 3579 PetscErrorCode ierr; 3580 3581 PetscFunctionBegin; 3582 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3583 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3584 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3585 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3586 PetscFunctionReturn(0); 3587 } 3588 3589 /*@C 3590 MatSeqAIJSetPreallocation - For good matrix assembly performance 3591 the user should preallocate the matrix storage by setting the parameter nz 3592 (or the array nnz). By setting these parameters accurately, performance 3593 during matrix assembly can be increased by more than a factor of 50. 3594 3595 Collective on MPI_Comm 3596 3597 Input Parameters: 3598 + B - The matrix 3599 . nz - number of nonzeros per row (same for all rows) 3600 - nnz - array containing the number of nonzeros in the various rows 3601 (possibly different for each row) or NULL 3602 3603 Notes: 3604 If nnz is given then nz is ignored 3605 3606 The AIJ format (also called the Yale sparse matrix format or 3607 compressed row storage), is fully compatible with standard Fortran 77 3608 storage. That is, the stored row and column indices can begin at 3609 either one (as in Fortran) or zero. See the users' manual for details. 3610 3611 Specify the preallocated storage with either nz or nnz (not both). 3612 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3613 allocation. For large problems you MUST preallocate memory or you 3614 will get TERRIBLE performance, see the users' manual chapter on matrices. 3615 3616 You can call MatGetInfo() to get information on how effective the preallocation was; 3617 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3618 You can also run with the option -info and look for messages with the string 3619 malloc in them to see if additional memory allocation was needed. 3620 3621 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3622 entries or columns indices 3623 3624 By default, this format uses inodes (identical nodes) when possible, to 3625 improve numerical efficiency of matrix-vector products and solves. We 3626 search for consecutive rows with the same nonzero structure, thereby 3627 reusing matrix information to achieve increased efficiency. 3628 3629 Options Database Keys: 3630 + -mat_no_inode - Do not use inodes 3631 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3632 3633 Level: intermediate 3634 3635 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3636 3637 @*/ 3638 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3639 { 3640 PetscErrorCode ierr; 3641 3642 PetscFunctionBegin; 3643 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3644 PetscValidType(B,1); 3645 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3646 PetscFunctionReturn(0); 3647 } 3648 3649 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3650 { 3651 Mat_SeqAIJ *b; 3652 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3653 PetscErrorCode ierr; 3654 PetscInt i; 3655 3656 PetscFunctionBegin; 3657 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3658 if (nz == MAT_SKIP_ALLOCATION) { 3659 skipallocation = PETSC_TRUE; 3660 nz = 0; 3661 } 3662 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3663 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3664 3665 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3666 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3667 if (nnz) { 3668 for (i=0; i<B->rmap->n; i++) { 3669 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]); 3670 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); 3671 } 3672 } 3673 3674 B->preallocated = PETSC_TRUE; 3675 3676 b = (Mat_SeqAIJ*)B->data; 3677 3678 if (!skipallocation) { 3679 if (!b->imax) { 3680 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3681 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3682 } 3683 if (!b->ipre) { 3684 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3685 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3686 } 3687 if (!nnz) { 3688 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3689 else if (nz < 0) nz = 1; 3690 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3691 nz = nz*B->rmap->n; 3692 } else { 3693 nz = 0; 3694 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3695 } 3696 /* b->ilen will count nonzeros in each row so far. */ 3697 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3698 3699 /* allocate the matrix space */ 3700 /* FIXME: should B's old memory be unlogged? */ 3701 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3702 if (B->structure_only) { 3703 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3704 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3705 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3706 } else { 3707 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3708 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3709 } 3710 b->i[0] = 0; 3711 for (i=1; i<B->rmap->n+1; i++) { 3712 b->i[i] = b->i[i-1] + b->imax[i-1]; 3713 } 3714 if (B->structure_only) { 3715 b->singlemalloc = PETSC_FALSE; 3716 b->free_a = PETSC_FALSE; 3717 } else { 3718 b->singlemalloc = PETSC_TRUE; 3719 b->free_a = PETSC_TRUE; 3720 } 3721 b->free_ij = PETSC_TRUE; 3722 } else { 3723 b->free_a = PETSC_FALSE; 3724 b->free_ij = PETSC_FALSE; 3725 } 3726 3727 if (b->ipre && nnz != b->ipre && b->imax) { 3728 /* reserve user-requested sparsity */ 3729 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3730 } 3731 3732 3733 b->nz = 0; 3734 b->maxnz = nz; 3735 B->info.nz_unneeded = (double)b->maxnz; 3736 if (realalloc) { 3737 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3738 } 3739 B->was_assembled = PETSC_FALSE; 3740 B->assembled = PETSC_FALSE; 3741 PetscFunctionReturn(0); 3742 } 3743 3744 3745 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3746 { 3747 Mat_SeqAIJ *a; 3748 PetscInt i; 3749 PetscErrorCode ierr; 3750 3751 PetscFunctionBegin; 3752 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3753 a = (Mat_SeqAIJ*)A->data; 3754 /* if no saved info, we error out */ 3755 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3756 3757 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"); 3758 3759 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3760 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3761 a->i[0] = 0; 3762 for (i=1; i<A->rmap->n+1; i++) { 3763 a->i[i] = a->i[i-1] + a->imax[i-1]; 3764 } 3765 A->preallocated = PETSC_TRUE; 3766 a->nz = 0; 3767 a->maxnz = a->i[A->rmap->n]; 3768 A->info.nz_unneeded = (double)a->maxnz; 3769 A->was_assembled = PETSC_FALSE; 3770 A->assembled = PETSC_FALSE; 3771 PetscFunctionReturn(0); 3772 } 3773 3774 /*@ 3775 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3776 3777 Input Parameters: 3778 + B - the matrix 3779 . i - the indices into j for the start of each row (starts with zero) 3780 . j - the column indices for each row (starts with zero) these must be sorted for each row 3781 - v - optional values in the matrix 3782 3783 Level: developer 3784 3785 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3786 3787 .keywords: matrix, aij, compressed row, sparse, sequential 3788 3789 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ 3790 @*/ 3791 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3792 { 3793 PetscErrorCode ierr; 3794 3795 PetscFunctionBegin; 3796 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3797 PetscValidType(B,1); 3798 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3799 PetscFunctionReturn(0); 3800 } 3801 3802 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3803 { 3804 PetscInt i; 3805 PetscInt m,n; 3806 PetscInt nz; 3807 PetscInt *nnz, nz_max = 0; 3808 PetscScalar *values; 3809 PetscErrorCode ierr; 3810 3811 PetscFunctionBegin; 3812 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3813 3814 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3815 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3816 3817 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3818 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3819 for (i = 0; i < m; i++) { 3820 nz = Ii[i+1]- Ii[i]; 3821 nz_max = PetscMax(nz_max, nz); 3822 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3823 nnz[i] = nz; 3824 } 3825 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3826 ierr = PetscFree(nnz);CHKERRQ(ierr); 3827 3828 if (v) { 3829 values = (PetscScalar*) v; 3830 } else { 3831 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3832 } 3833 3834 for (i = 0; i < m; i++) { 3835 nz = Ii[i+1] - Ii[i]; 3836 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3837 } 3838 3839 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3840 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3841 3842 if (!v) { 3843 ierr = PetscFree(values);CHKERRQ(ierr); 3844 } 3845 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3846 PetscFunctionReturn(0); 3847 } 3848 3849 #include <../src/mat/impls/dense/seq/dense.h> 3850 #include <petsc/private/kernels/petscaxpy.h> 3851 3852 /* 3853 Computes (B'*A')' since computing B*A directly is untenable 3854 3855 n p p 3856 ( ) ( ) ( ) 3857 m ( A ) * n ( B ) = m ( C ) 3858 ( ) ( ) ( ) 3859 3860 */ 3861 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3862 { 3863 PetscErrorCode ierr; 3864 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3865 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3866 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3867 PetscInt i,n,m,q,p; 3868 const PetscInt *ii,*idx; 3869 const PetscScalar *b,*a,*a_q; 3870 PetscScalar *c,*c_q; 3871 3872 PetscFunctionBegin; 3873 m = A->rmap->n; 3874 n = A->cmap->n; 3875 p = B->cmap->n; 3876 a = sub_a->v; 3877 b = sub_b->a; 3878 c = sub_c->v; 3879 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3880 3881 ii = sub_b->i; 3882 idx = sub_b->j; 3883 for (i=0; i<n; i++) { 3884 q = ii[i+1] - ii[i]; 3885 while (q-->0) { 3886 c_q = c + m*(*idx); 3887 a_q = a + m*i; 3888 PetscKernelAXPY(c_q,*b,a_q,m); 3889 idx++; 3890 b++; 3891 } 3892 } 3893 PetscFunctionReturn(0); 3894 } 3895 3896 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3897 { 3898 PetscErrorCode ierr; 3899 PetscInt m=A->rmap->n,n=B->cmap->n; 3900 Mat Cmat; 3901 3902 PetscFunctionBegin; 3903 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); 3904 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3905 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3906 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3907 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3908 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3909 3910 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3911 3912 *C = Cmat; 3913 PetscFunctionReturn(0); 3914 } 3915 3916 /* ----------------------------------------------------------------*/ 3917 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3918 { 3919 PetscErrorCode ierr; 3920 3921 PetscFunctionBegin; 3922 if (scall == MAT_INITIAL_MATRIX) { 3923 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3924 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3925 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3926 } 3927 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3928 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3929 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3930 PetscFunctionReturn(0); 3931 } 3932 3933 3934 /*MC 3935 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3936 based on compressed sparse row format. 3937 3938 Options Database Keys: 3939 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3940 3941 Level: beginner 3942 3943 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3944 M*/ 3945 3946 /*MC 3947 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3948 3949 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3950 and MATMPIAIJ otherwise. As a result, for single process communicators, 3951 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3952 for communicators controlling multiple processes. It is recommended that you call both of 3953 the above preallocation routines for simplicity. 3954 3955 Options Database Keys: 3956 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3957 3958 Developer Notes: 3959 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 3960 enough exist. 3961 3962 Level: beginner 3963 3964 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3965 M*/ 3966 3967 /*MC 3968 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3969 3970 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3971 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3972 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3973 for communicators controlling multiple processes. It is recommended that you call both of 3974 the above preallocation routines for simplicity. 3975 3976 Options Database Keys: 3977 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3978 3979 Level: beginner 3980 3981 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3982 M*/ 3983 3984 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3985 #if defined(PETSC_HAVE_ELEMENTAL) 3986 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3987 #endif 3988 #if defined(PETSC_HAVE_HYPRE) 3989 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3990 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3991 #endif 3992 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3993 3994 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3995 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3996 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3997 #endif 3998 3999 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4000 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4001 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 4002 4003 /*@C 4004 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 4005 4006 Not Collective 4007 4008 Input Parameter: 4009 . mat - a MATSEQAIJ matrix 4010 4011 Output Parameter: 4012 . array - pointer to the data 4013 4014 Level: intermediate 4015 4016 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4017 @*/ 4018 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4019 { 4020 PetscErrorCode ierr; 4021 4022 PetscFunctionBegin; 4023 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4024 PetscFunctionReturn(0); 4025 } 4026 4027 /*@C 4028 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4029 4030 Not Collective 4031 4032 Input Parameter: 4033 . mat - a MATSEQAIJ matrix 4034 4035 Output Parameter: 4036 . nz - the maximum number of nonzeros in any row 4037 4038 Level: intermediate 4039 4040 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4041 @*/ 4042 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4043 { 4044 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4045 4046 PetscFunctionBegin; 4047 *nz = aij->rmax; 4048 PetscFunctionReturn(0); 4049 } 4050 4051 /*@C 4052 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4053 4054 Not Collective 4055 4056 Input Parameters: 4057 . mat - a MATSEQAIJ matrix 4058 . array - pointer to the data 4059 4060 Level: intermediate 4061 4062 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4063 @*/ 4064 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4065 { 4066 PetscErrorCode ierr; 4067 4068 PetscFunctionBegin; 4069 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4070 PetscFunctionReturn(0); 4071 } 4072 4073 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4074 { 4075 Mat_SeqAIJ *b; 4076 PetscErrorCode ierr; 4077 PetscMPIInt size; 4078 4079 PetscFunctionBegin; 4080 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4081 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4082 4083 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4084 4085 B->data = (void*)b; 4086 4087 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4088 4089 b->row = 0; 4090 b->col = 0; 4091 b->icol = 0; 4092 b->reallocs = 0; 4093 b->ignorezeroentries = PETSC_FALSE; 4094 b->roworiented = PETSC_TRUE; 4095 b->nonew = 0; 4096 b->diag = 0; 4097 b->solve_work = 0; 4098 B->spptr = 0; 4099 b->saved_values = 0; 4100 b->idiag = 0; 4101 b->mdiag = 0; 4102 b->ssor_work = 0; 4103 b->omega = 1.0; 4104 b->fshift = 0.0; 4105 b->idiagvalid = PETSC_FALSE; 4106 b->ibdiagvalid = PETSC_FALSE; 4107 b->keepnonzeropattern = PETSC_FALSE; 4108 4109 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4110 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4111 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4112 4113 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4114 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4115 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4116 #endif 4117 4118 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4119 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4120 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4121 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4122 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4123 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4124 #if defined(PETSC_HAVE_MKL_SPARSE) 4125 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4126 #endif 4127 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4128 #if defined(PETSC_HAVE_ELEMENTAL) 4129 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4130 #endif 4131 #if defined(PETSC_HAVE_HYPRE) 4132 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4133 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4134 #endif 4135 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4136 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4137 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4138 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4139 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4140 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4141 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4142 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4143 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4144 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4145 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4146 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4147 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaij_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 4148 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4149 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4150 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4151 PetscFunctionReturn(0); 4152 } 4153 4154 /* 4155 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4156 */ 4157 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4158 { 4159 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4160 PetscErrorCode ierr; 4161 PetscInt i,m = A->rmap->n; 4162 4163 PetscFunctionBegin; 4164 c = (Mat_SeqAIJ*)C->data; 4165 4166 C->factortype = A->factortype; 4167 c->row = 0; 4168 c->col = 0; 4169 c->icol = 0; 4170 c->reallocs = 0; 4171 4172 C->assembled = PETSC_TRUE; 4173 4174 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4175 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4176 4177 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4178 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4179 for (i=0; i<m; i++) { 4180 c->imax[i] = a->imax[i]; 4181 c->ilen[i] = a->ilen[i]; 4182 } 4183 4184 /* allocate the matrix space */ 4185 if (mallocmatspace) { 4186 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4187 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4188 4189 c->singlemalloc = PETSC_TRUE; 4190 4191 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4192 if (m > 0) { 4193 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4194 if (cpvalues == MAT_COPY_VALUES) { 4195 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4196 } else { 4197 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4198 } 4199 } 4200 } 4201 4202 c->ignorezeroentries = a->ignorezeroentries; 4203 c->roworiented = a->roworiented; 4204 c->nonew = a->nonew; 4205 if (a->diag) { 4206 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4207 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4208 for (i=0; i<m; i++) { 4209 c->diag[i] = a->diag[i]; 4210 } 4211 } else c->diag = 0; 4212 4213 c->solve_work = 0; 4214 c->saved_values = 0; 4215 c->idiag = 0; 4216 c->ssor_work = 0; 4217 c->keepnonzeropattern = a->keepnonzeropattern; 4218 c->free_a = PETSC_TRUE; 4219 c->free_ij = PETSC_TRUE; 4220 4221 c->rmax = a->rmax; 4222 c->nz = a->nz; 4223 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4224 C->preallocated = PETSC_TRUE; 4225 4226 c->compressedrow.use = a->compressedrow.use; 4227 c->compressedrow.nrows = a->compressedrow.nrows; 4228 if (a->compressedrow.use) { 4229 i = a->compressedrow.nrows; 4230 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4231 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4232 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4233 } else { 4234 c->compressedrow.use = PETSC_FALSE; 4235 c->compressedrow.i = NULL; 4236 c->compressedrow.rindex = NULL; 4237 } 4238 c->nonzerorowcnt = a->nonzerorowcnt; 4239 C->nonzerostate = A->nonzerostate; 4240 4241 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4242 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4243 PetscFunctionReturn(0); 4244 } 4245 4246 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4247 { 4248 PetscErrorCode ierr; 4249 4250 PetscFunctionBegin; 4251 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4252 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4253 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4254 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4255 } 4256 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4257 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4258 PetscFunctionReturn(0); 4259 } 4260 4261 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4262 { 4263 Mat_SeqAIJ *a; 4264 PetscErrorCode ierr; 4265 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4266 int fd; 4267 PetscMPIInt size; 4268 MPI_Comm comm; 4269 PetscInt bs = newMat->rmap->bs; 4270 4271 PetscFunctionBegin; 4272 /* force binary viewer to load .info file if it has not yet done so */ 4273 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4274 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4275 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4276 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4277 4278 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4279 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4280 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4281 if (bs < 0) bs = 1; 4282 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4283 4284 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4285 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4286 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4287 M = header[1]; N = header[2]; nz = header[3]; 4288 4289 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4290 4291 /* read in row lengths */ 4292 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4293 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4294 4295 /* check if sum of rowlengths is same as nz */ 4296 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4297 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); 4298 4299 /* set global size if not set already*/ 4300 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4301 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4302 } else { 4303 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4304 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4305 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4306 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4307 } 4308 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); 4309 } 4310 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4311 a = (Mat_SeqAIJ*)newMat->data; 4312 4313 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4314 4315 /* read in nonzero values */ 4316 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4317 4318 /* set matrix "i" values */ 4319 a->i[0] = 0; 4320 for (i=1; i<= M; i++) { 4321 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4322 a->ilen[i-1] = rowlengths[i-1]; 4323 } 4324 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4325 4326 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4327 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4328 PetscFunctionReturn(0); 4329 } 4330 4331 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4332 { 4333 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4334 PetscErrorCode ierr; 4335 #if defined(PETSC_USE_COMPLEX) 4336 PetscInt k; 4337 #endif 4338 4339 PetscFunctionBegin; 4340 /* If the matrix dimensions are not equal,or no of nonzeros */ 4341 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4342 *flg = PETSC_FALSE; 4343 PetscFunctionReturn(0); 4344 } 4345 4346 /* if the a->i are the same */ 4347 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4348 if (!*flg) PetscFunctionReturn(0); 4349 4350 /* if a->j are the same */ 4351 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4352 if (!*flg) PetscFunctionReturn(0); 4353 4354 /* if a->a are the same */ 4355 #if defined(PETSC_USE_COMPLEX) 4356 for (k=0; k<a->nz; k++) { 4357 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4358 *flg = PETSC_FALSE; 4359 PetscFunctionReturn(0); 4360 } 4361 } 4362 #else 4363 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4364 #endif 4365 PetscFunctionReturn(0); 4366 } 4367 4368 /*@ 4369 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4370 provided by the user. 4371 4372 Collective on MPI_Comm 4373 4374 Input Parameters: 4375 + comm - must be an MPI communicator of size 1 4376 . m - number of rows 4377 . n - number of columns 4378 . i - row indices 4379 . j - column indices 4380 - a - matrix values 4381 4382 Output Parameter: 4383 . mat - the matrix 4384 4385 Level: intermediate 4386 4387 Notes: 4388 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4389 once the matrix is destroyed and not before 4390 4391 You cannot set new nonzero locations into this matrix, that will generate an error. 4392 4393 The i and j indices are 0 based 4394 4395 The format which is used for the sparse matrix input, is equivalent to a 4396 row-major ordering.. i.e for the following matrix, the input data expected is 4397 as shown 4398 4399 $ 1 0 0 4400 $ 2 0 3 4401 $ 4 5 6 4402 $ 4403 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4404 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4405 $ v = {1,2,3,4,5,6} [size = 6] 4406 4407 4408 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4409 4410 @*/ 4411 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4412 { 4413 PetscErrorCode ierr; 4414 PetscInt ii; 4415 Mat_SeqAIJ *aij; 4416 #if defined(PETSC_USE_DEBUG) 4417 PetscInt jj; 4418 #endif 4419 4420 PetscFunctionBegin; 4421 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4422 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4423 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4424 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4425 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4426 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4427 aij = (Mat_SeqAIJ*)(*mat)->data; 4428 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4429 4430 aij->i = i; 4431 aij->j = j; 4432 aij->a = a; 4433 aij->singlemalloc = PETSC_FALSE; 4434 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4435 aij->free_a = PETSC_FALSE; 4436 aij->free_ij = PETSC_FALSE; 4437 4438 for (ii=0; ii<m; ii++) { 4439 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4440 #if defined(PETSC_USE_DEBUG) 4441 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]); 4442 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4443 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); 4444 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); 4445 } 4446 #endif 4447 } 4448 #if defined(PETSC_USE_DEBUG) 4449 for (ii=0; ii<aij->i[m]; ii++) { 4450 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4451 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]); 4452 } 4453 #endif 4454 4455 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4456 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4457 PetscFunctionReturn(0); 4458 } 4459 /*@C 4460 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4461 provided by the user. 4462 4463 Collective on MPI_Comm 4464 4465 Input Parameters: 4466 + comm - must be an MPI communicator of size 1 4467 . m - number of rows 4468 . n - number of columns 4469 . i - row indices 4470 . j - column indices 4471 . a - matrix values 4472 . nz - number of nonzeros 4473 - idx - 0 or 1 based 4474 4475 Output Parameter: 4476 . mat - the matrix 4477 4478 Level: intermediate 4479 4480 Notes: 4481 The i and j indices are 0 based 4482 4483 The format which is used for the sparse matrix input, is equivalent to a 4484 row-major ordering.. i.e for the following matrix, the input data expected is 4485 as shown: 4486 4487 1 0 0 4488 2 0 3 4489 4 5 6 4490 4491 i = {0,1,1,2,2,2} 4492 j = {0,0,2,0,1,2} 4493 v = {1,2,3,4,5,6} 4494 4495 4496 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4497 4498 @*/ 4499 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4500 { 4501 PetscErrorCode ierr; 4502 PetscInt ii, *nnz, one = 1,row,col; 4503 4504 4505 PetscFunctionBegin; 4506 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4507 for (ii = 0; ii < nz; ii++) { 4508 nnz[i[ii] - !!idx] += 1; 4509 } 4510 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4511 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4512 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4513 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4514 for (ii = 0; ii < nz; ii++) { 4515 if (idx) { 4516 row = i[ii] - 1; 4517 col = j[ii] - 1; 4518 } else { 4519 row = i[ii]; 4520 col = j[ii]; 4521 } 4522 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4523 } 4524 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4525 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4526 ierr = PetscFree(nnz);CHKERRQ(ierr); 4527 PetscFunctionReturn(0); 4528 } 4529 4530 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4531 { 4532 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4533 PetscErrorCode ierr; 4534 4535 PetscFunctionBegin; 4536 a->idiagvalid = PETSC_FALSE; 4537 a->ibdiagvalid = PETSC_FALSE; 4538 4539 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4540 PetscFunctionReturn(0); 4541 } 4542 4543 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4544 { 4545 PetscErrorCode ierr; 4546 PetscMPIInt size; 4547 4548 PetscFunctionBegin; 4549 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4550 if (size == 1) { 4551 if (scall == MAT_INITIAL_MATRIX) { 4552 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4553 } else { 4554 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4555 } 4556 } else { 4557 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4558 } 4559 PetscFunctionReturn(0); 4560 } 4561 4562 /* 4563 Permute A into C's *local* index space using rowemb,colemb. 4564 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4565 of [0,m), colemb is in [0,n). 4566 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4567 */ 4568 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4569 { 4570 /* If making this function public, change the error returned in this function away from _PLIB. */ 4571 PetscErrorCode ierr; 4572 Mat_SeqAIJ *Baij; 4573 PetscBool seqaij; 4574 PetscInt m,n,*nz,i,j,count; 4575 PetscScalar v; 4576 const PetscInt *rowindices,*colindices; 4577 4578 PetscFunctionBegin; 4579 if (!B) PetscFunctionReturn(0); 4580 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4581 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4582 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4583 if (rowemb) { 4584 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4585 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); 4586 } else { 4587 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4588 } 4589 if (colemb) { 4590 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4591 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); 4592 } else { 4593 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4594 } 4595 4596 Baij = (Mat_SeqAIJ*)(B->data); 4597 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4598 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4599 for (i=0; i<B->rmap->n; i++) { 4600 nz[i] = Baij->i[i+1] - Baij->i[i]; 4601 } 4602 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4603 ierr = PetscFree(nz);CHKERRQ(ierr); 4604 } 4605 if (pattern == SUBSET_NONZERO_PATTERN) { 4606 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4607 } 4608 count = 0; 4609 rowindices = NULL; 4610 colindices = NULL; 4611 if (rowemb) { 4612 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4613 } 4614 if (colemb) { 4615 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4616 } 4617 for (i=0; i<B->rmap->n; i++) { 4618 PetscInt row; 4619 row = i; 4620 if (rowindices) row = rowindices[i]; 4621 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4622 PetscInt col; 4623 col = Baij->j[count]; 4624 if (colindices) col = colindices[col]; 4625 v = Baij->a[count]; 4626 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4627 ++count; 4628 } 4629 } 4630 /* FIXME: set C's nonzerostate correctly. */ 4631 /* Assembly for C is necessary. */ 4632 C->preallocated = PETSC_TRUE; 4633 C->assembled = PETSC_TRUE; 4634 C->was_assembled = PETSC_FALSE; 4635 PetscFunctionReturn(0); 4636 } 4637 4638 PetscFunctionList MatSeqAIJList = NULL; 4639 4640 /*@C 4641 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4642 4643 Collective on Mat 4644 4645 Input Parameters: 4646 + mat - the matrix object 4647 - matype - matrix type 4648 4649 Options Database Key: 4650 . -mat_seqai_type <method> - for example seqaijcrl 4651 4652 4653 Level: intermediate 4654 4655 .keywords: Mat, MatType, set, method 4656 4657 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4658 @*/ 4659 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4660 { 4661 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4662 PetscBool sametype; 4663 4664 PetscFunctionBegin; 4665 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4666 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4667 if (sametype) PetscFunctionReturn(0); 4668 4669 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4670 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4671 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4672 PetscFunctionReturn(0); 4673 } 4674 4675 4676 /*@C 4677 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4678 4679 Not Collective 4680 4681 Input Parameters: 4682 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4683 - function - routine to convert to subtype 4684 4685 Notes: 4686 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4687 4688 4689 Then, your matrix can be chosen with the procedural interface at runtime via the option 4690 $ -mat_seqaij_type my_mat 4691 4692 Level: advanced 4693 4694 .keywords: Mat, register 4695 4696 .seealso: MatSeqAIJRegisterAll() 4697 4698 4699 Level: advanced 4700 @*/ 4701 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4702 { 4703 PetscErrorCode ierr; 4704 4705 PetscFunctionBegin; 4706 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4707 PetscFunctionReturn(0); 4708 } 4709 4710 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4711 4712 /*@C 4713 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4714 4715 Not Collective 4716 4717 Level: advanced 4718 4719 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4720 4721 .keywords: KSP, register, all 4722 4723 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4724 @*/ 4725 PetscErrorCode MatSeqAIJRegisterAll(void) 4726 { 4727 PetscErrorCode ierr; 4728 4729 PetscFunctionBegin; 4730 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4731 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4732 4733 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4734 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4735 #if defined(PETSC_HAVE_MKL_SPARSE) 4736 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4737 #endif 4738 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4739 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4740 #endif 4741 PetscFunctionReturn(0); 4742 } 4743 4744 /* 4745 Special version for direct calls from Fortran 4746 */ 4747 #include <petsc/private/fortranimpl.h> 4748 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4749 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4750 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4751 #define matsetvaluesseqaij_ matsetvaluesseqaij 4752 #endif 4753 4754 /* Change these macros so can be used in void function */ 4755 #undef CHKERRQ 4756 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4757 #undef SETERRQ2 4758 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4759 #undef SETERRQ3 4760 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4761 4762 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) 4763 { 4764 Mat A = *AA; 4765 PetscInt m = *mm, n = *nn; 4766 InsertMode is = *isis; 4767 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4768 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4769 PetscInt *imax,*ai,*ailen; 4770 PetscErrorCode ierr; 4771 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4772 MatScalar *ap,value,*aa; 4773 PetscBool ignorezeroentries = a->ignorezeroentries; 4774 PetscBool roworiented = a->roworiented; 4775 4776 PetscFunctionBegin; 4777 MatCheckPreallocated(A,1); 4778 imax = a->imax; 4779 ai = a->i; 4780 ailen = a->ilen; 4781 aj = a->j; 4782 aa = a->a; 4783 4784 for (k=0; k<m; k++) { /* loop over added rows */ 4785 row = im[k]; 4786 if (row < 0) continue; 4787 #if defined(PETSC_USE_DEBUG) 4788 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4789 #endif 4790 rp = aj + ai[row]; ap = aa + ai[row]; 4791 rmax = imax[row]; nrow = ailen[row]; 4792 low = 0; 4793 high = nrow; 4794 for (l=0; l<n; l++) { /* loop over added columns */ 4795 if (in[l] < 0) continue; 4796 #if defined(PETSC_USE_DEBUG) 4797 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4798 #endif 4799 col = in[l]; 4800 if (roworiented) value = v[l + k*n]; 4801 else value = v[k + l*m]; 4802 4803 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4804 4805 if (col <= lastcol) low = 0; 4806 else high = nrow; 4807 lastcol = col; 4808 while (high-low > 5) { 4809 t = (low+high)/2; 4810 if (rp[t] > col) high = t; 4811 else low = t; 4812 } 4813 for (i=low; i<high; i++) { 4814 if (rp[i] > col) break; 4815 if (rp[i] == col) { 4816 if (is == ADD_VALUES) ap[i] += value; 4817 else ap[i] = value; 4818 goto noinsert; 4819 } 4820 } 4821 if (value == 0.0 && ignorezeroentries) goto noinsert; 4822 if (nonew == 1) goto noinsert; 4823 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4824 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4825 N = nrow++ - 1; a->nz++; high++; 4826 /* shift up all the later entries in this row */ 4827 for (ii=N; ii>=i; ii--) { 4828 rp[ii+1] = rp[ii]; 4829 ap[ii+1] = ap[ii]; 4830 } 4831 rp[i] = col; 4832 ap[i] = value; 4833 A->nonzerostate++; 4834 noinsert:; 4835 low = i + 1; 4836 } 4837 ailen[row] = nrow; 4838 } 4839 PetscFunctionReturnVoid(); 4840 } 4841