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