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