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 #if defined(PETSC_THREADCOMM_ACTIVE) 13 #include <petscthreadcomm.h> 14 #endif 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 18 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 19 { 20 PetscErrorCode ierr; 21 PetscInt i,m,n; 22 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 23 24 PetscFunctionBegin; 25 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 26 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 27 if (type == NORM_2) { 28 for (i=0; i<aij->i[m]; i++) { 29 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 30 } 31 } else if (type == NORM_1) { 32 for (i=0; i<aij->i[m]; i++) { 33 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 34 } 35 } else if (type == NORM_INFINITY) { 36 for (i=0; i<aij->i[m]; i++) { 37 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 38 } 39 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 40 41 if (type == NORM_2) { 42 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatFindOffBlockDiagonalEntries_SeqAIJ" 49 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 50 { 51 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 52 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 53 const PetscInt *jj = a->j,*ii = a->i; 54 PetscInt *rows; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 for (i=0; i<m; i++) { 59 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 60 cnt++; 61 } 62 } 63 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 64 cnt = 0; 65 for (i=0; i<m; i++) { 66 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 67 rows[cnt] = i; 68 cnt++; 69 } 70 } 71 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ_Private" 77 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 78 { 79 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 80 const MatScalar *aa = a->a; 81 PetscInt i,m=A->rmap->n,cnt = 0; 82 const PetscInt *jj = a->j,*diag; 83 PetscInt *rows; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 88 diag = a->diag; 89 for (i=0; i<m; i++) { 90 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 91 cnt++; 92 } 93 } 94 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 95 cnt = 0; 96 for (i=0; i<m; i++) { 97 if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 98 rows[cnt++] = i; 99 } 100 } 101 *nrows = cnt; 102 *zrows = rows; 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ" 108 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 109 { 110 PetscInt nrows,*rows; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 *zrows = NULL; 115 ierr = MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);CHKERRQ(ierr); 116 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatFindNonzeroRows_SeqAIJ" 122 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 123 { 124 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 125 const MatScalar *aa; 126 PetscInt m=A->rmap->n,cnt = 0; 127 const PetscInt *ii; 128 PetscInt n,i,j,*rows; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 *keptrows = 0; 133 ii = a->i; 134 for (i=0; i<m; i++) { 135 n = ii[i+1] - ii[i]; 136 if (!n) { 137 cnt++; 138 goto ok1; 139 } 140 aa = a->a + ii[i]; 141 for (j=0; j<n; j++) { 142 if (aa[j] != 0.0) goto ok1; 143 } 144 cnt++; 145 ok1:; 146 } 147 if (!cnt) PetscFunctionReturn(0); 148 ierr = PetscMalloc1((A->rmap->n-cnt),&rows);CHKERRQ(ierr); 149 cnt = 0; 150 for (i=0; i<m; i++) { 151 n = ii[i+1] - ii[i]; 152 if (!n) continue; 153 aa = a->a + ii[i]; 154 for (j=0; j<n; j++) { 155 if (aa[j] != 0.0) { 156 rows[cnt++] = i; 157 break; 158 } 159 } 160 } 161 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 167 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 168 { 169 PetscErrorCode ierr; 170 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 171 PetscInt i,*diag, m = Y->rmap->n; 172 MatScalar *aa = aij->a; 173 PetscScalar *v; 174 PetscBool missing; 175 176 PetscFunctionBegin; 177 if (Y->assembled) { 178 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);CHKERRQ(ierr); 179 if (!missing) { 180 diag = aij->diag; 181 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 182 if (is == INSERT_VALUES) { 183 for (i=0; i<m; i++) { 184 aa[diag[i]] = v[i]; 185 } 186 } else { 187 for (i=0; i<m; i++) { 188 aa[diag[i]] += v[i]; 189 } 190 } 191 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 195 } 196 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 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,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 #undef __FUNCT__ 235 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 236 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 237 { 238 PetscErrorCode ierr; 239 240 PetscFunctionBegin; 241 if (!ia) PetscFunctionReturn(0); 242 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 243 ierr = PetscFree(*ia);CHKERRQ(ierr); 244 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 245 } 246 PetscFunctionReturn(0); 247 } 248 249 #undef __FUNCT__ 250 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 251 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 252 { 253 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 254 PetscErrorCode ierr; 255 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 256 PetscInt nz = a->i[m],row,*jj,mr,col; 257 258 PetscFunctionBegin; 259 *nn = n; 260 if (!ia) PetscFunctionReturn(0); 261 if (symmetric) { 262 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 263 } else { 264 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 265 ierr = PetscMalloc1((n+1),&cia);CHKERRQ(ierr); 266 ierr = PetscMalloc1((nz+1),&cja);CHKERRQ(ierr); 267 jj = a->j; 268 for (i=0; i<nz; i++) { 269 collengths[jj[i]]++; 270 } 271 cia[0] = oshift; 272 for (i=0; i<n; i++) { 273 cia[i+1] = cia[i] + collengths[i]; 274 } 275 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 276 jj = a->j; 277 for (row=0; row<m; row++) { 278 mr = a->i[row+1] - a->i[row]; 279 for (i=0; i<mr; i++) { 280 col = *jj++; 281 282 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 283 } 284 } 285 ierr = PetscFree(collengths);CHKERRQ(ierr); 286 *ia = cia; *ja = cja; 287 } 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 293 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 294 { 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 if (!ia) PetscFunctionReturn(0); 299 300 ierr = PetscFree(*ia);CHKERRQ(ierr); 301 ierr = PetscFree(*ja);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 /* 306 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 307 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 308 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 309 */ 310 #undef __FUNCT__ 311 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 312 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 313 { 314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 315 PetscErrorCode ierr; 316 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 317 PetscInt nz = a->i[m],row,*jj,mr,col; 318 PetscInt *cspidx; 319 320 PetscFunctionBegin; 321 *nn = n; 322 if (!ia) PetscFunctionReturn(0); 323 324 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 325 ierr = PetscMalloc1((n+1),&cia);CHKERRQ(ierr); 326 ierr = PetscMalloc1((nz+1),&cja);CHKERRQ(ierr); 327 ierr = PetscMalloc1((nz+1),&cspidx);CHKERRQ(ierr); 328 jj = a->j; 329 for (i=0; i<nz; i++) { 330 collengths[jj[i]]++; 331 } 332 cia[0] = oshift; 333 for (i=0; i<n; i++) { 334 cia[i+1] = cia[i] + collengths[i]; 335 } 336 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 337 jj = a->j; 338 for (row=0; row<m; row++) { 339 mr = a->i[row+1] - a->i[row]; 340 for (i=0; i<mr; i++) { 341 col = *jj++; 342 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 343 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 344 } 345 } 346 ierr = PetscFree(collengths);CHKERRQ(ierr); 347 *ia = cia; *ja = cja; 348 *spidx = cspidx; 349 PetscFunctionReturn(0); 350 } 351 352 #undef __FUNCT__ 353 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 354 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 355 { 356 PetscErrorCode ierr; 357 358 PetscFunctionBegin; 359 ierr = MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);CHKERRQ(ierr); 360 ierr = PetscFree(*spidx);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 366 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 367 { 368 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 369 PetscInt *ai = a->i; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 374 PetscFunctionReturn(0); 375 } 376 377 /* 378 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 379 380 - a single row of values is set with each call 381 - no row or column indices are negative or (in error) larger than the number of rows or columns 382 - the values are always added to the matrix, not set 383 - no new locations are introduced in the nonzero structure of the matrix 384 385 This does NOT assume the global column indices are sorted 386 387 */ 388 389 #include <petsc-private/isimpl.h> 390 #undef __FUNCT__ 391 #define __FUNCT__ "MatSeqAIJSetValuesLocalFast" 392 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 393 { 394 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 395 PetscInt low,high,t,row,nrow,i,col,l; 396 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 397 PetscInt lastcol = -1; 398 MatScalar *ap,value,*aa = a->a; 399 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 400 401 row = ridx[im[0]]; 402 rp = aj + ai[row]; 403 ap = aa + ai[row]; 404 nrow = ailen[row]; 405 low = 0; 406 high = nrow; 407 for (l=0; l<n; l++) { /* loop over added columns */ 408 col = cidx[in[l]]; 409 value = v[l]; 410 411 if (col <= lastcol) low = 0; 412 else high = nrow; 413 lastcol = col; 414 while (high-low > 5) { 415 t = (low+high)/2; 416 if (rp[t] > col) high = t; 417 else low = t; 418 } 419 for (i=low; i<high; i++) { 420 if (rp[i] == col) { 421 ap[i] += value; 422 low = i + 1; 423 break; 424 } 425 } 426 } 427 return 0; 428 } 429 430 #undef __FUNCT__ 431 #define __FUNCT__ "MatSetValues_SeqAIJ" 432 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 433 { 434 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 435 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 436 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 437 PetscErrorCode ierr; 438 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 439 MatScalar *ap,value,*aa = a->a; 440 PetscBool ignorezeroentries = a->ignorezeroentries; 441 PetscBool roworiented = a->roworiented; 442 443 PetscFunctionBegin; 444 if (v) PetscValidScalarPointer(v,6); 445 for (k=0; k<m; k++) { /* loop over added rows */ 446 row = im[k]; 447 if (row < 0) continue; 448 #if defined(PETSC_USE_DEBUG) 449 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); 450 #endif 451 rp = aj + ai[row]; ap = aa + ai[row]; 452 rmax = imax[row]; nrow = ailen[row]; 453 low = 0; 454 high = nrow; 455 for (l=0; l<n; l++) { /* loop over added columns */ 456 if (in[l] < 0) continue; 457 #if defined(PETSC_USE_DEBUG) 458 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); 459 #endif 460 col = in[l]; 461 if (v) { 462 if (roworiented) { 463 value = v[l + k*n]; 464 } else { 465 value = v[k + l*m]; 466 } 467 } else { 468 value = 0.; 469 } 470 if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES)) continue; 471 472 if (col <= lastcol) low = 0; 473 else high = nrow; 474 lastcol = col; 475 while (high-low > 5) { 476 t = (low+high)/2; 477 if (rp[t] > col) high = t; 478 else low = t; 479 } 480 for (i=low; i<high; i++) { 481 if (rp[i] > col) break; 482 if (rp[i] == col) { 483 if (is == ADD_VALUES) ap[i] += value; 484 else ap[i] = value; 485 low = i + 1; 486 goto noinsert; 487 } 488 } 489 if (value == 0.0 && ignorezeroentries) goto noinsert; 490 if (nonew == 1) goto noinsert; 491 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 492 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 493 N = nrow++ - 1; a->nz++; high++; 494 /* shift up all the later entries in this row */ 495 for (ii=N; ii>=i; ii--) { 496 rp[ii+1] = rp[ii]; 497 ap[ii+1] = ap[ii]; 498 } 499 rp[i] = col; 500 ap[i] = value; 501 low = i + 1; 502 A->nonzerostate++; 503 noinsert:; 504 } 505 ailen[row] = nrow; 506 } 507 PetscFunctionReturn(0); 508 } 509 510 511 #undef __FUNCT__ 512 #define __FUNCT__ "MatGetValues_SeqAIJ" 513 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 514 { 515 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 516 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 517 PetscInt *ai = a->i,*ailen = a->ilen; 518 MatScalar *ap,*aa = a->a; 519 520 PetscFunctionBegin; 521 for (k=0; k<m; k++) { /* loop over rows */ 522 row = im[k]; 523 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 524 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); 525 rp = aj + ai[row]; ap = aa + ai[row]; 526 nrow = ailen[row]; 527 for (l=0; l<n; l++) { /* loop over columns */ 528 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 529 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); 530 col = in[l]; 531 high = nrow; low = 0; /* assume unsorted */ 532 while (high-low > 5) { 533 t = (low+high)/2; 534 if (rp[t] > col) high = t; 535 else low = t; 536 } 537 for (i=low; i<high; i++) { 538 if (rp[i] > col) break; 539 if (rp[i] == col) { 540 *v++ = ap[i]; 541 goto finished; 542 } 543 } 544 *v++ = 0.0; 545 finished:; 546 } 547 } 548 PetscFunctionReturn(0); 549 } 550 551 552 #undef __FUNCT__ 553 #define __FUNCT__ "MatView_SeqAIJ_Binary" 554 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 555 { 556 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 557 PetscErrorCode ierr; 558 PetscInt i,*col_lens; 559 int fd; 560 FILE *file; 561 562 PetscFunctionBegin; 563 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 564 ierr = PetscMalloc1((4+A->rmap->n),&col_lens);CHKERRQ(ierr); 565 566 col_lens[0] = MAT_FILE_CLASSID; 567 col_lens[1] = A->rmap->n; 568 col_lens[2] = A->cmap->n; 569 col_lens[3] = a->nz; 570 571 /* store lengths of each row and write (including header) to file */ 572 for (i=0; i<A->rmap->n; i++) { 573 col_lens[4+i] = a->i[i+1] - a->i[i]; 574 } 575 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 576 ierr = PetscFree(col_lens);CHKERRQ(ierr); 577 578 /* store column indices (zero start index) */ 579 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 580 581 /* store nonzero values */ 582 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 583 584 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 585 if (file) { 586 fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs)); 587 } 588 PetscFunctionReturn(0); 589 } 590 591 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 592 593 #undef __FUNCT__ 594 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 595 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 596 { 597 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 598 PetscErrorCode ierr; 599 PetscInt i,j,m = A->rmap->n; 600 const char *name; 601 PetscViewerFormat format; 602 603 PetscFunctionBegin; 604 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 605 if (format == PETSC_VIEWER_ASCII_MATLAB) { 606 PetscInt nofinalvalue = 0; 607 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 608 /* Need a dummy value to ensure the dimension of the matrix. */ 609 nofinalvalue = 1; 610 } 611 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 612 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 613 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 614 #if defined(PETSC_USE_COMPLEX) 615 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,4);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 616 #else 617 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 618 #endif 619 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 620 621 for (i=0; i<m; i++) { 622 for (j=a->i[i]; j<a->i[i+1]; j++) { 623 #if defined(PETSC_USE_COMPLEX) 624 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); 625 #else 626 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);CHKERRQ(ierr); 627 #endif 628 } 629 } 630 if (nofinalvalue) { 631 #if defined(PETSC_USE_COMPLEX) 632 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",m,A->cmap->n,0.,0.);CHKERRQ(ierr); 633 #else 634 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 635 #endif 636 } 637 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 638 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 639 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 640 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 641 PetscFunctionReturn(0); 642 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 643 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 644 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 645 for (i=0; i<m; i++) { 646 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 647 for (j=a->i[i]; j<a->i[i+1]; j++) { 648 #if defined(PETSC_USE_COMPLEX) 649 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 650 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 651 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 652 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 653 } else if (PetscRealPart(a->a[j]) != 0.0) { 654 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 655 } 656 #else 657 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr);} 658 #endif 659 } 660 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 661 } 662 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 663 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 664 PetscInt nzd=0,fshift=1,*sptr; 665 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 666 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 667 ierr = PetscMalloc1((m+1),&sptr);CHKERRQ(ierr); 668 for (i=0; i<m; i++) { 669 sptr[i] = nzd+1; 670 for (j=a->i[i]; j<a->i[i+1]; j++) { 671 if (a->j[j] >= i) { 672 #if defined(PETSC_USE_COMPLEX) 673 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 674 #else 675 if (a->a[j] != 0.0) nzd++; 676 #endif 677 } 678 } 679 } 680 sptr[m] = nzd+1; 681 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 682 for (i=0; i<m+1; i+=6) { 683 if (i+4<m) { 684 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); 685 } else if (i+3<m) { 686 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); 687 } else if (i+2<m) { 688 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr); 689 } else if (i+1<m) { 690 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr); 691 } else if (i<m) { 692 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr); 693 } else { 694 ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr); 695 } 696 } 697 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 698 ierr = PetscFree(sptr);CHKERRQ(ierr); 699 for (i=0; i<m; i++) { 700 for (j=a->i[i]; j<a->i[i+1]; j++) { 701 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 702 } 703 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 704 } 705 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 706 for (i=0; i<m; i++) { 707 for (j=a->i[i]; j<a->i[i+1]; j++) { 708 if (a->j[j] >= i) { 709 #if defined(PETSC_USE_COMPLEX) 710 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 711 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 712 } 713 #else 714 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);CHKERRQ(ierr);} 715 #endif 716 } 717 } 718 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 719 } 720 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 721 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 722 PetscInt cnt = 0,jcnt; 723 PetscScalar value; 724 #if defined(PETSC_USE_COMPLEX) 725 PetscBool realonly = PETSC_TRUE; 726 727 for (i=0; i<a->i[m]; i++) { 728 if (PetscImaginaryPart(a->a[i]) != 0.0) { 729 realonly = PETSC_FALSE; 730 break; 731 } 732 } 733 #endif 734 735 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 736 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 737 for (i=0; i<m; i++) { 738 jcnt = 0; 739 for (j=0; j<A->cmap->n; j++) { 740 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 741 value = a->a[cnt++]; 742 jcnt++; 743 } else { 744 value = 0.0; 745 } 746 #if defined(PETSC_USE_COMPLEX) 747 if (realonly) { 748 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));CHKERRQ(ierr); 749 } else { 750 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));CHKERRQ(ierr); 751 } 752 #else 753 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);CHKERRQ(ierr); 754 #endif 755 } 756 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 757 } 758 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 759 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 760 PetscInt fshift=1; 761 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 762 #if defined(PETSC_USE_COMPLEX) 763 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 764 #else 765 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 766 #endif 767 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 768 for (i=0; i<m; i++) { 769 for (j=a->i[i]; j<a->i[i+1]; j++) { 770 #if defined(PETSC_USE_COMPLEX) 771 if (PetscImaginaryPart(a->a[j]) > 0.0) { 772 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); 773 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 774 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); 775 } else { 776 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 777 } 778 #else 779 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);CHKERRQ(ierr); 780 #endif 781 } 782 } 783 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 784 } else { 785 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 786 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 787 if (A->factortype) { 788 for (i=0; i<m; i++) { 789 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 790 /* L part */ 791 for (j=a->i[i]; j<a->i[i+1]; j++) { 792 #if defined(PETSC_USE_COMPLEX) 793 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 if (PetscImaginaryPart(a->a[j]) < 0.0) { 796 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 797 } else { 798 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 799 } 800 #else 801 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 802 #endif 803 } 804 /* diagonal */ 805 j = a->diag[i]; 806 #if defined(PETSC_USE_COMPLEX) 807 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 if (PetscImaginaryPart(a->a[j]) < 0.0) { 810 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); 811 } else { 812 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 813 } 814 #else 815 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)(1.0/a->a[j]));CHKERRQ(ierr); 816 #endif 817 818 /* U part */ 819 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 820 #if defined(PETSC_USE_COMPLEX) 821 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 if (PetscImaginaryPart(a->a[j]) < 0.0) { 824 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 825 } else { 826 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 827 } 828 #else 829 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 830 #endif 831 } 832 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 833 } 834 } else { 835 for (i=0; i<m; i++) { 836 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 837 for (j=a->i[i]; j<a->i[i+1]; j++) { 838 #if defined(PETSC_USE_COMPLEX) 839 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 if (PetscImaginaryPart(a->a[j]) < 0.0) { 842 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 843 } else { 844 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 845 } 846 #else 847 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 848 #endif 849 } 850 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 851 } 852 } 853 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 854 } 855 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 856 PetscFunctionReturn(0); 857 } 858 859 #include <petscdraw.h> 860 #undef __FUNCT__ 861 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 862 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 863 { 864 Mat A = (Mat) Aa; 865 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 866 PetscErrorCode ierr; 867 PetscInt i,j,m = A->rmap->n,color; 868 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 869 PetscViewer viewer; 870 PetscViewerFormat format; 871 872 PetscFunctionBegin; 873 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 874 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 875 876 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 877 /* loop over matrix elements drawing boxes */ 878 879 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 880 /* Blue for negative, Cyan for zero and Red for positive */ 881 color = PETSC_DRAW_BLUE; 882 for (i=0; i<m; i++) { 883 y_l = m - i - 1.0; y_r = y_l + 1.0; 884 for (j=a->i[i]; j<a->i[i+1]; j++) { 885 x_l = a->j[j]; x_r = x_l + 1.0; 886 if (PetscRealPart(a->a[j]) >= 0.) continue; 887 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 888 } 889 } 890 color = PETSC_DRAW_CYAN; 891 for (i=0; i<m; i++) { 892 y_l = m - i - 1.0; y_r = y_l + 1.0; 893 for (j=a->i[i]; j<a->i[i+1]; j++) { 894 x_l = a->j[j]; x_r = x_l + 1.0; 895 if (a->a[j] != 0.) continue; 896 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 897 } 898 } 899 color = PETSC_DRAW_RED; 900 for (i=0; i<m; i++) { 901 y_l = m - i - 1.0; y_r = y_l + 1.0; 902 for (j=a->i[i]; j<a->i[i+1]; j++) { 903 x_l = a->j[j]; x_r = x_l + 1.0; 904 if (PetscRealPart(a->a[j]) <= 0.) continue; 905 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 906 } 907 } 908 } else { 909 /* use contour shading to indicate magnitude of values */ 910 /* first determine max of all nonzero values */ 911 PetscInt nz = a->nz,count; 912 PetscDraw popup; 913 PetscReal scale; 914 915 for (i=0; i<nz; i++) { 916 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 917 } 918 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 919 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 920 if (popup) { 921 ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr); 922 } 923 count = 0; 924 for (i=0; i<m; i++) { 925 y_l = m - i - 1.0; y_r = y_l + 1.0; 926 for (j=a->i[i]; j<a->i[i+1]; j++) { 927 x_l = a->j[j]; x_r = x_l + 1.0; 928 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 929 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 930 count++; 931 } 932 } 933 } 934 PetscFunctionReturn(0); 935 } 936 937 #include <petscdraw.h> 938 #undef __FUNCT__ 939 #define __FUNCT__ "MatView_SeqAIJ_Draw" 940 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 941 { 942 PetscErrorCode ierr; 943 PetscDraw draw; 944 PetscReal xr,yr,xl,yl,h,w; 945 PetscBool isnull; 946 947 PetscFunctionBegin; 948 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 949 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 950 if (isnull) PetscFunctionReturn(0); 951 952 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 953 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 954 xr += w; yr += h; xl = -w; yl = -h; 955 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 956 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 957 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 958 PetscFunctionReturn(0); 959 } 960 961 #undef __FUNCT__ 962 #define __FUNCT__ "MatView_SeqAIJ" 963 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 964 { 965 PetscErrorCode ierr; 966 PetscBool iascii,isbinary,isdraw; 967 968 PetscFunctionBegin; 969 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 970 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 971 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 972 if (iascii) { 973 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 974 } else if (isbinary) { 975 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 976 } else if (isdraw) { 977 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 978 } 979 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 980 PetscFunctionReturn(0); 981 } 982 983 #undef __FUNCT__ 984 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 985 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 986 { 987 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 988 PetscErrorCode ierr; 989 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 990 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 991 MatScalar *aa = a->a,*ap; 992 PetscReal ratio = 0.6; 993 994 PetscFunctionBegin; 995 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 996 997 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 998 for (i=1; i<m; i++) { 999 /* move each row back by the amount of empty slots (fshift) before it*/ 1000 fshift += imax[i-1] - ailen[i-1]; 1001 rmax = PetscMax(rmax,ailen[i]); 1002 if (fshift) { 1003 ip = aj + ai[i]; 1004 ap = aa + ai[i]; 1005 N = ailen[i]; 1006 for (j=0; j<N; j++) { 1007 ip[j-fshift] = ip[j]; 1008 ap[j-fshift] = ap[j]; 1009 } 1010 } 1011 ai[i] = ai[i-1] + ailen[i-1]; 1012 } 1013 if (m) { 1014 fshift += imax[m-1] - ailen[m-1]; 1015 ai[m] = ai[m-1] + ailen[m-1]; 1016 } 1017 1018 /* reset ilen and imax for each row */ 1019 a->nonzerorowcnt = 0; 1020 for (i=0; i<m; i++) { 1021 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1022 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 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 ierr = MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 1038 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 1039 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1040 PetscFunctionReturn(0); 1041 } 1042 1043 #undef __FUNCT__ 1044 #define __FUNCT__ "MatRealPart_SeqAIJ" 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 #undef __FUNCT__ 1059 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 1060 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1061 { 1062 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1063 PetscInt i,nz = a->nz; 1064 MatScalar *aa = a->a; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1069 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #if defined(PETSC_THREADCOMM_ACTIVE) 1074 PetscErrorCode MatZeroEntries_SeqAIJ_Kernel(PetscInt thread_id,Mat A) 1075 { 1076 PetscErrorCode ierr; 1077 PetscInt *trstarts=A->rmap->trstarts; 1078 PetscInt n,start,end; 1079 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1080 1081 start = trstarts[thread_id]; 1082 end = trstarts[thread_id+1]; 1083 n = a->i[end] - a->i[start]; 1084 ierr = PetscMemzero(a->a+a->i[start],n*sizeof(PetscScalar));CHKERRQ(ierr); 1085 return 0; 1086 } 1087 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1090 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1091 { 1092 PetscErrorCode ierr; 1093 1094 PetscFunctionBegin; 1095 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatZeroEntries_SeqAIJ_Kernel,1,A);CHKERRQ(ierr); 1096 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1097 PetscFunctionReturn(0); 1098 } 1099 #else 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1102 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1103 { 1104 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1105 PetscErrorCode ierr; 1106 1107 PetscFunctionBegin; 1108 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1109 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 #endif 1113 1114 #undef __FUNCT__ 1115 #define __FUNCT__ "MatDestroy_SeqAIJ" 1116 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1117 { 1118 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1119 PetscErrorCode ierr; 1120 1121 PetscFunctionBegin; 1122 #if defined(PETSC_USE_LOG) 1123 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1124 #endif 1125 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1126 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1127 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1128 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1129 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1130 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1131 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1132 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1133 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1134 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1135 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1136 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 1137 ierr = MatDestroy(&a->XtoY);CHKERRQ(ierr); 1138 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1139 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1140 1141 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1142 ierr = PetscFree(A->data);CHKERRQ(ierr); 1143 1144 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1145 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1146 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1147 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1148 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1149 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1150 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1151 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1152 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1153 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1154 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 #undef __FUNCT__ 1159 #define __FUNCT__ "MatSetOption_SeqAIJ" 1160 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1161 { 1162 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1163 PetscErrorCode ierr; 1164 1165 PetscFunctionBegin; 1166 switch (op) { 1167 case MAT_ROW_ORIENTED: 1168 a->roworiented = flg; 1169 break; 1170 case MAT_KEEP_NONZERO_PATTERN: 1171 a->keepnonzeropattern = flg; 1172 break; 1173 case MAT_NEW_NONZERO_LOCATIONS: 1174 a->nonew = (flg ? 0 : 1); 1175 break; 1176 case MAT_NEW_NONZERO_LOCATION_ERR: 1177 a->nonew = (flg ? -1 : 0); 1178 break; 1179 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1180 a->nonew = (flg ? -2 : 0); 1181 break; 1182 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1183 a->nounused = (flg ? -1 : 0); 1184 break; 1185 case MAT_IGNORE_ZERO_ENTRIES: 1186 a->ignorezeroentries = flg; 1187 break; 1188 case MAT_SPD: 1189 case MAT_SYMMETRIC: 1190 case MAT_STRUCTURALLY_SYMMETRIC: 1191 case MAT_HERMITIAN: 1192 case MAT_SYMMETRY_ETERNAL: 1193 /* These options are handled directly by MatSetOption() */ 1194 break; 1195 case MAT_NEW_DIAGONALS: 1196 case MAT_IGNORE_OFF_PROC_ENTRIES: 1197 case MAT_USE_HASH_TABLE: 1198 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1199 break; 1200 case MAT_USE_INODES: 1201 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1202 break; 1203 default: 1204 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1205 } 1206 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 1212 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1213 { 1214 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1215 PetscErrorCode ierr; 1216 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1217 PetscScalar *aa=a->a,*x,zero=0.0; 1218 1219 PetscFunctionBegin; 1220 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1221 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1222 1223 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1224 PetscInt *diag=a->diag; 1225 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1226 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1227 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 ierr = VecSet(v,zero);CHKERRQ(ierr); 1232 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1233 for (i=0; i<n; i++) { 1234 nz = ai[i+1] - ai[i]; 1235 if (!nz) x[i] = 0.0; 1236 for (j=ai[i]; j<ai[i+1]; j++) { 1237 if (aj[j] == i) { 1238 x[i] = aa[j]; 1239 break; 1240 } 1241 } 1242 } 1243 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1244 PetscFunctionReturn(0); 1245 } 1246 1247 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1248 #undef __FUNCT__ 1249 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 1250 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1251 { 1252 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1253 PetscScalar *x,*y; 1254 PetscErrorCode ierr; 1255 PetscInt m = A->rmap->n; 1256 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1257 MatScalar *v; 1258 PetscScalar alpha; 1259 PetscInt n,i,j,*idx,*ii,*ridx=NULL; 1260 Mat_CompressedRow cprow = a->compressedrow; 1261 PetscBool usecprow = cprow.use; 1262 #endif 1263 1264 PetscFunctionBegin; 1265 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1266 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1267 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1268 1269 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1270 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1271 #else 1272 if (usecprow) { 1273 m = cprow.nrows; 1274 ii = cprow.i; 1275 ridx = cprow.rindex; 1276 } else { 1277 ii = a->i; 1278 } 1279 for (i=0; i<m; i++) { 1280 idx = a->j + ii[i]; 1281 v = a->a + ii[i]; 1282 n = ii[i+1] - ii[i]; 1283 if (usecprow) { 1284 alpha = x[ridx[i]]; 1285 } else { 1286 alpha = x[i]; 1287 } 1288 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1289 } 1290 #endif 1291 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1292 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1293 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1294 PetscFunctionReturn(0); 1295 } 1296 1297 #undef __FUNCT__ 1298 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 1299 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1300 { 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1305 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1306 PetscFunctionReturn(0); 1307 } 1308 1309 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1310 #if defined(PETSC_THREADCOMM_ACTIVE) 1311 PetscErrorCode MatMult_SeqAIJ_Kernel(PetscInt thread_id,Mat A,Vec xx,Vec yy) 1312 { 1313 PetscErrorCode ierr; 1314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1315 PetscScalar *y; 1316 const PetscScalar *x; 1317 const MatScalar *aa; 1318 PetscInt *trstarts=A->rmap->trstarts; 1319 PetscInt n,start,end,i; 1320 const PetscInt *aj,*ai; 1321 PetscScalar sum; 1322 1323 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1324 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1325 start = trstarts[thread_id]; 1326 end = trstarts[thread_id+1]; 1327 aj = a->j; 1328 aa = a->a; 1329 ai = a->i; 1330 for (i=start; i<end; i++) { 1331 n = ai[i+1] - ai[i]; 1332 aj = a->j + ai[i]; 1333 aa = a->a + ai[i]; 1334 sum = 0.0; 1335 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1336 y[i] = sum; 1337 } 1338 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1339 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1340 return 0; 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "MatMult_SeqAIJ" 1345 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1346 { 1347 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1348 PetscScalar *y; 1349 const PetscScalar *x; 1350 const MatScalar *aa; 1351 PetscErrorCode ierr; 1352 PetscInt m=A->rmap->n; 1353 const PetscInt *aj,*ii,*ridx=NULL; 1354 PetscInt n,i; 1355 PetscScalar sum; 1356 PetscBool usecprow=a->compressedrow.use; 1357 1358 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1359 #pragma disjoint(*x,*y,*aa) 1360 #endif 1361 1362 PetscFunctionBegin; 1363 aj = a->j; 1364 aa = a->a; 1365 ii = a->i; 1366 if (usecprow) { /* use compressed row format */ 1367 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1368 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1369 m = a->compressedrow.nrows; 1370 ii = a->compressedrow.i; 1371 ridx = a->compressedrow.rindex; 1372 for (i=0; i<m; i++) { 1373 n = ii[i+1] - ii[i]; 1374 aj = a->j + ii[i]; 1375 aa = a->a + ii[i]; 1376 sum = 0.0; 1377 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1378 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1379 y[*ridx++] = sum; 1380 } 1381 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1382 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1383 } else { /* do not use compressed row format */ 1384 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1385 fortranmultaij_(&m,x,ii,aj,aa,y); 1386 #else 1387 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1388 #endif 1389 } 1390 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1391 PetscFunctionReturn(0); 1392 } 1393 #else 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "MatMult_SeqAIJ" 1396 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1397 { 1398 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1399 PetscScalar *y; 1400 const PetscScalar *x; 1401 const MatScalar *aa; 1402 PetscErrorCode ierr; 1403 PetscInt m=A->rmap->n; 1404 const PetscInt *aj,*ii,*ridx=NULL; 1405 PetscInt n,i; 1406 PetscScalar sum; 1407 PetscBool usecprow=a->compressedrow.use; 1408 1409 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1410 #pragma disjoint(*x,*y,*aa) 1411 #endif 1412 1413 PetscFunctionBegin; 1414 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1415 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1416 aj = a->j; 1417 aa = a->a; 1418 ii = a->i; 1419 if (usecprow) { /* use compressed row format */ 1420 m = a->compressedrow.nrows; 1421 ii = a->compressedrow.i; 1422 ridx = a->compressedrow.rindex; 1423 for (i=0; i<m; i++) { 1424 n = ii[i+1] - ii[i]; 1425 aj = a->j + ii[i]; 1426 aa = a->a + ii[i]; 1427 sum = 0.0; 1428 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1429 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1430 y[*ridx++] = sum; 1431 } 1432 } else { /* do not use compressed row format */ 1433 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1434 fortranmultaij_(&m,x,ii,aj,aa,y); 1435 #else 1436 #if defined(PETSC_THREADCOMM_ACTIVE) 1437 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1438 #else 1439 for (i=0; i<m; i++) { 1440 n = ii[i+1] - ii[i]; 1441 aj = a->j + ii[i]; 1442 aa = a->a + ii[i]; 1443 sum = 0.0; 1444 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1445 y[i] = sum; 1446 } 1447 #endif 1448 #endif 1449 } 1450 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1451 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1452 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 #endif 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatMultMax_SeqAIJ" 1459 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1460 { 1461 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1462 PetscScalar *y; 1463 const PetscScalar *x; 1464 const MatScalar *aa; 1465 PetscErrorCode ierr; 1466 PetscInt m=A->rmap->n; 1467 const PetscInt *aj,*ii,*ridx=NULL; 1468 PetscInt n,i,nonzerorow=0; 1469 PetscScalar sum; 1470 PetscBool usecprow=a->compressedrow.use; 1471 1472 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1473 #pragma disjoint(*x,*y,*aa) 1474 #endif 1475 1476 PetscFunctionBegin; 1477 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1478 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1479 aj = a->j; 1480 aa = a->a; 1481 ii = a->i; 1482 if (usecprow) { /* use compressed row format */ 1483 m = a->compressedrow.nrows; 1484 ii = a->compressedrow.i; 1485 ridx = a->compressedrow.rindex; 1486 for (i=0; i<m; i++) { 1487 n = ii[i+1] - ii[i]; 1488 aj = a->j + ii[i]; 1489 aa = a->a + ii[i]; 1490 sum = 0.0; 1491 nonzerorow += (n>0); 1492 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1493 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1494 y[*ridx++] = sum; 1495 } 1496 } else { /* do not use compressed row format */ 1497 for (i=0; i<m; i++) { 1498 n = ii[i+1] - ii[i]; 1499 aj = a->j + ii[i]; 1500 aa = a->a + ii[i]; 1501 sum = 0.0; 1502 nonzerorow += (n>0); 1503 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1504 y[i] = sum; 1505 } 1506 } 1507 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1508 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1509 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 #undef __FUNCT__ 1514 #define __FUNCT__ "MatMultAddMax_SeqAIJ" 1515 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1516 { 1517 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1518 PetscScalar *y,*z; 1519 const PetscScalar *x; 1520 const MatScalar *aa; 1521 PetscErrorCode ierr; 1522 PetscInt m = A->rmap->n,*aj,*ii; 1523 PetscInt n,i,*ridx=NULL; 1524 PetscScalar sum; 1525 PetscBool usecprow=a->compressedrow.use; 1526 1527 PetscFunctionBegin; 1528 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1529 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1530 if (zz != yy) { 1531 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1532 } else { 1533 z = y; 1534 } 1535 1536 aj = a->j; 1537 aa = a->a; 1538 ii = a->i; 1539 if (usecprow) { /* use compressed row format */ 1540 if (zz != yy) { 1541 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1542 } 1543 m = a->compressedrow.nrows; 1544 ii = a->compressedrow.i; 1545 ridx = a->compressedrow.rindex; 1546 for (i=0; i<m; i++) { 1547 n = ii[i+1] - ii[i]; 1548 aj = a->j + ii[i]; 1549 aa = a->a + ii[i]; 1550 sum = y[*ridx]; 1551 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1552 z[*ridx++] = sum; 1553 } 1554 } else { /* do not use compressed row format */ 1555 for (i=0; i<m; i++) { 1556 n = ii[i+1] - ii[i]; 1557 aj = a->j + ii[i]; 1558 aa = a->a + ii[i]; 1559 sum = y[i]; 1560 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1561 z[i] = sum; 1562 } 1563 } 1564 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1565 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1566 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1567 if (zz != yy) { 1568 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1569 } 1570 PetscFunctionReturn(0); 1571 } 1572 1573 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1576 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1577 { 1578 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1579 PetscScalar *y,*z; 1580 const PetscScalar *x; 1581 const MatScalar *aa; 1582 PetscErrorCode ierr; 1583 PetscInt m = A->rmap->n,*aj,*ii; 1584 PetscInt n,i,*ridx=NULL; 1585 PetscScalar sum; 1586 PetscBool usecprow=a->compressedrow.use; 1587 1588 PetscFunctionBegin; 1589 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1590 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1591 if (zz != yy) { 1592 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1593 } else { 1594 z = y; 1595 } 1596 1597 aj = a->j; 1598 aa = a->a; 1599 ii = a->i; 1600 if (usecprow) { /* use compressed row format */ 1601 if (zz != yy) { 1602 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1603 } 1604 m = a->compressedrow.nrows; 1605 ii = a->compressedrow.i; 1606 ridx = a->compressedrow.rindex; 1607 for (i=0; i<m; i++) { 1608 n = ii[i+1] - ii[i]; 1609 aj = a->j + ii[i]; 1610 aa = a->a + ii[i]; 1611 sum = y[*ridx]; 1612 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1613 z[*ridx++] = sum; 1614 } 1615 } else { /* do not use compressed row format */ 1616 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1617 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1618 #else 1619 for (i=0; i<m; i++) { 1620 n = ii[i+1] - ii[i]; 1621 aj = a->j + ii[i]; 1622 aa = a->a + ii[i]; 1623 sum = y[i]; 1624 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1625 z[i] = sum; 1626 } 1627 #endif 1628 } 1629 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1630 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1631 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1632 if (zz != yy) { 1633 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1634 } 1635 #if defined(PETSC_HAVE_CUSP) 1636 /* 1637 ierr = VecView(xx,0);CHKERRQ(ierr); 1638 ierr = VecView(zz,0);CHKERRQ(ierr); 1639 ierr = MatView(A,0);CHKERRQ(ierr); 1640 */ 1641 #endif 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /* 1646 Adds diagonal pointers to sparse matrix structure. 1647 */ 1648 #undef __FUNCT__ 1649 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1650 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1651 { 1652 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1653 PetscErrorCode ierr; 1654 PetscInt i,j,m = A->rmap->n; 1655 1656 PetscFunctionBegin; 1657 if (!a->diag) { 1658 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1659 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1660 } 1661 for (i=0; i<A->rmap->n; i++) { 1662 a->diag[i] = a->i[i+1]; 1663 for (j=a->i[i]; j<a->i[i+1]; j++) { 1664 if (a->j[j] == i) { 1665 a->diag[i] = j; 1666 break; 1667 } 1668 } 1669 } 1670 PetscFunctionReturn(0); 1671 } 1672 1673 /* 1674 Checks for missing diagonals 1675 */ 1676 #undef __FUNCT__ 1677 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1678 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1679 { 1680 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1681 PetscInt *diag,*jj = a->j,i; 1682 1683 PetscFunctionBegin; 1684 *missing = PETSC_FALSE; 1685 if (A->rmap->n > 0 && !jj) { 1686 *missing = PETSC_TRUE; 1687 if (d) *d = 0; 1688 PetscInfo(A,"Matrix has no entries therefore is missing diagonal"); 1689 } else { 1690 diag = a->diag; 1691 for (i=0; i<A->rmap->n; i++) { 1692 if (jj[diag[i]] != i) { 1693 *missing = PETSC_TRUE; 1694 if (d) *d = i; 1695 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1696 break; 1697 } 1698 } 1699 } 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1705 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1706 { 1707 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1708 PetscErrorCode ierr; 1709 PetscInt i,*diag,m = A->rmap->n; 1710 MatScalar *v = a->a; 1711 PetscScalar *idiag,*mdiag; 1712 1713 PetscFunctionBegin; 1714 if (a->idiagvalid) PetscFunctionReturn(0); 1715 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1716 diag = a->diag; 1717 if (!a->idiag) { 1718 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1719 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1720 v = a->a; 1721 } 1722 mdiag = a->mdiag; 1723 idiag = a->idiag; 1724 1725 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1726 for (i=0; i<m; i++) { 1727 mdiag[i] = v[diag[i]]; 1728 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1729 idiag[i] = 1.0/v[diag[i]]; 1730 } 1731 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1732 } else { 1733 for (i=0; i<m; i++) { 1734 mdiag[i] = v[diag[i]]; 1735 idiag[i] = omega/(fshift + v[diag[i]]); 1736 } 1737 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1738 } 1739 a->idiagvalid = PETSC_TRUE; 1740 PetscFunctionReturn(0); 1741 } 1742 1743 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1744 #undef __FUNCT__ 1745 #define __FUNCT__ "MatSOR_SeqAIJ" 1746 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1747 { 1748 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1749 PetscScalar *x,d,sum,*t,scale; 1750 const MatScalar *v = a->a,*idiag=0,*mdiag; 1751 const PetscScalar *b, *bs,*xb, *ts; 1752 PetscErrorCode ierr; 1753 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1754 const PetscInt *idx,*diag; 1755 1756 PetscFunctionBegin; 1757 its = its*lits; 1758 1759 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1760 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1761 a->fshift = fshift; 1762 a->omega = omega; 1763 1764 diag = a->diag; 1765 t = a->ssor_work; 1766 idiag = a->idiag; 1767 mdiag = a->mdiag; 1768 1769 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1770 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1771 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1772 if (flag == SOR_APPLY_UPPER) { 1773 /* apply (U + D/omega) to the vector */ 1774 bs = b; 1775 for (i=0; i<m; i++) { 1776 d = fshift + mdiag[i]; 1777 n = a->i[i+1] - diag[i] - 1; 1778 idx = a->j + diag[i] + 1; 1779 v = a->a + diag[i] + 1; 1780 sum = b[i]*d/omega; 1781 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1782 x[i] = sum; 1783 } 1784 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1785 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1786 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1787 PetscFunctionReturn(0); 1788 } 1789 1790 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1791 else if (flag & SOR_EISENSTAT) { 1792 /* Let A = L + U + D; where L is lower trianglar, 1793 U is upper triangular, E = D/omega; This routine applies 1794 1795 (L + E)^{-1} A (U + E)^{-1} 1796 1797 to a vector efficiently using Eisenstat's trick. 1798 */ 1799 scale = (2.0/omega) - 1.0; 1800 1801 /* x = (E + U)^{-1} b */ 1802 for (i=m-1; i>=0; i--) { 1803 n = a->i[i+1] - diag[i] - 1; 1804 idx = a->j + diag[i] + 1; 1805 v = a->a + diag[i] + 1; 1806 sum = b[i]; 1807 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1808 x[i] = sum*idiag[i]; 1809 } 1810 1811 /* t = b - (2*E - D)x */ 1812 v = a->a; 1813 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1814 1815 /* t = (E + L)^{-1}t */ 1816 ts = t; 1817 diag = a->diag; 1818 for (i=0; i<m; i++) { 1819 n = diag[i] - a->i[i]; 1820 idx = a->j + a->i[i]; 1821 v = a->a + a->i[i]; 1822 sum = t[i]; 1823 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1824 t[i] = sum*idiag[i]; 1825 /* x = x + t */ 1826 x[i] += t[i]; 1827 } 1828 1829 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1830 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1831 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1832 PetscFunctionReturn(0); 1833 } 1834 if (flag & SOR_ZERO_INITIAL_GUESS) { 1835 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1836 for (i=0; i<m; i++) { 1837 n = diag[i] - a->i[i]; 1838 idx = a->j + a->i[i]; 1839 v = a->a + a->i[i]; 1840 sum = b[i]; 1841 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1842 t[i] = sum; 1843 x[i] = sum*idiag[i]; 1844 } 1845 xb = t; 1846 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1847 } else xb = b; 1848 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1849 for (i=m-1; i>=0; i--) { 1850 n = a->i[i+1] - diag[i] - 1; 1851 idx = a->j + diag[i] + 1; 1852 v = a->a + diag[i] + 1; 1853 sum = xb[i]; 1854 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1855 if (xb == b) { 1856 x[i] = sum*idiag[i]; 1857 } else { 1858 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1859 } 1860 } 1861 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1862 } 1863 its--; 1864 } 1865 while (its--) { 1866 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1867 for (i=0; i<m; i++) { 1868 /* lower */ 1869 n = diag[i] - a->i[i]; 1870 idx = a->j + a->i[i]; 1871 v = a->a + a->i[i]; 1872 sum = b[i]; 1873 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1874 t[i] = sum; /* save application of the lower-triangular part */ 1875 /* upper */ 1876 n = a->i[i+1] - diag[i] - 1; 1877 idx = a->j + diag[i] + 1; 1878 v = a->a + diag[i] + 1; 1879 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1880 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1881 } 1882 xb = t; 1883 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1884 } else xb = b; 1885 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1886 for (i=m-1; i>=0; i--) { 1887 sum = xb[i]; 1888 if (xb == b) { 1889 /* whole matrix (no checkpointing available) */ 1890 n = a->i[i+1] - a->i[i]; 1891 idx = a->j + a->i[i]; 1892 v = a->a + a->i[i]; 1893 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1894 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1895 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1896 n = a->i[i+1] - diag[i] - 1; 1897 idx = a->j + diag[i] + 1; 1898 v = a->a + diag[i] + 1; 1899 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1900 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1901 } 1902 } 1903 if (xb == b) { 1904 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1905 } else { 1906 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1907 } 1908 } 1909 } 1910 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1911 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 1916 #undef __FUNCT__ 1917 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1918 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1919 { 1920 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1921 1922 PetscFunctionBegin; 1923 info->block_size = 1.0; 1924 info->nz_allocated = (double)a->maxnz; 1925 info->nz_used = (double)a->nz; 1926 info->nz_unneeded = (double)(a->maxnz - a->nz); 1927 info->assemblies = (double)A->num_ass; 1928 info->mallocs = (double)A->info.mallocs; 1929 info->memory = ((PetscObject)A)->mem; 1930 if (A->factortype) { 1931 info->fill_ratio_given = A->info.fill_ratio_given; 1932 info->fill_ratio_needed = A->info.fill_ratio_needed; 1933 info->factor_mallocs = A->info.factor_mallocs; 1934 } else { 1935 info->fill_ratio_given = 0; 1936 info->fill_ratio_needed = 0; 1937 info->factor_mallocs = 0; 1938 } 1939 PetscFunctionReturn(0); 1940 } 1941 1942 #undef __FUNCT__ 1943 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1944 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1945 { 1946 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1947 PetscInt i,m = A->rmap->n - 1,d = 0; 1948 PetscErrorCode ierr; 1949 const PetscScalar *xx; 1950 PetscScalar *bb; 1951 PetscBool missing; 1952 1953 PetscFunctionBegin; 1954 if (x && b) { 1955 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1956 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1957 for (i=0; i<N; i++) { 1958 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1959 bb[rows[i]] = diag*xx[rows[i]]; 1960 } 1961 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1962 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1963 } 1964 1965 if (a->keepnonzeropattern) { 1966 for (i=0; i<N; i++) { 1967 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1968 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1969 } 1970 if (diag != 0.0) { 1971 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1972 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1973 for (i=0; i<N; i++) { 1974 a->a[a->diag[rows[i]]] = diag; 1975 } 1976 } 1977 } else { 1978 if (diag != 0.0) { 1979 for (i=0; i<N; i++) { 1980 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1981 if (a->ilen[rows[i]] > 0) { 1982 a->ilen[rows[i]] = 1; 1983 a->a[a->i[rows[i]]] = diag; 1984 a->j[a->i[rows[i]]] = rows[i]; 1985 } else { /* in case row was completely empty */ 1986 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1987 } 1988 } 1989 } else { 1990 for (i=0; i<N; i++) { 1991 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1992 a->ilen[rows[i]] = 0; 1993 } 1994 } 1995 A->nonzerostate++; 1996 } 1997 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1998 PetscFunctionReturn(0); 1999 } 2000 2001 #undef __FUNCT__ 2002 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 2003 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2004 { 2005 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2006 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2007 PetscErrorCode ierr; 2008 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2009 const PetscScalar *xx; 2010 PetscScalar *bb; 2011 2012 PetscFunctionBegin; 2013 if (x && b) { 2014 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2015 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2016 vecs = PETSC_TRUE; 2017 } 2018 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2019 for (i=0; i<N; i++) { 2020 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2021 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2022 2023 zeroed[rows[i]] = PETSC_TRUE; 2024 } 2025 for (i=0; i<A->rmap->n; i++) { 2026 if (!zeroed[i]) { 2027 for (j=a->i[i]; j<a->i[i+1]; j++) { 2028 if (zeroed[a->j[j]]) { 2029 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2030 a->a[j] = 0.0; 2031 } 2032 } 2033 } else if (vecs) bb[i] = diag*xx[i]; 2034 } 2035 if (x && b) { 2036 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2037 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2038 } 2039 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2040 if (diag != 0.0) { 2041 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2042 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 2043 for (i=0; i<N; i++) { 2044 a->a[a->diag[rows[i]]] = diag; 2045 } 2046 } 2047 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "MatGetRow_SeqAIJ" 2053 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2054 { 2055 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2056 PetscInt *itmp; 2057 2058 PetscFunctionBegin; 2059 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2060 2061 *nz = a->i[row+1] - a->i[row]; 2062 if (v) *v = a->a + a->i[row]; 2063 if (idx) { 2064 itmp = a->j + a->i[row]; 2065 if (*nz) *idx = itmp; 2066 else *idx = 0; 2067 } 2068 PetscFunctionReturn(0); 2069 } 2070 2071 /* remove this function? */ 2072 #undef __FUNCT__ 2073 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 2074 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2075 { 2076 PetscFunctionBegin; 2077 PetscFunctionReturn(0); 2078 } 2079 2080 #undef __FUNCT__ 2081 #define __FUNCT__ "MatNorm_SeqAIJ" 2082 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2083 { 2084 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2085 MatScalar *v = a->a; 2086 PetscReal sum = 0.0; 2087 PetscErrorCode ierr; 2088 PetscInt i,j; 2089 2090 PetscFunctionBegin; 2091 if (type == NORM_FROBENIUS) { 2092 for (i=0; i<a->nz; i++) { 2093 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2094 } 2095 *nrm = PetscSqrtReal(sum); 2096 } else if (type == NORM_1) { 2097 PetscReal *tmp; 2098 PetscInt *jj = a->j; 2099 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2100 *nrm = 0.0; 2101 for (j=0; j<a->nz; j++) { 2102 tmp[*jj++] += PetscAbsScalar(*v); v++; 2103 } 2104 for (j=0; j<A->cmap->n; j++) { 2105 if (tmp[j] > *nrm) *nrm = tmp[j]; 2106 } 2107 ierr = PetscFree(tmp);CHKERRQ(ierr); 2108 } else if (type == NORM_INFINITY) { 2109 *nrm = 0.0; 2110 for (j=0; j<A->rmap->n; j++) { 2111 v = a->a + a->i[j]; 2112 sum = 0.0; 2113 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2114 sum += PetscAbsScalar(*v); v++; 2115 } 2116 if (sum > *nrm) *nrm = sum; 2117 } 2118 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2123 #undef __FUNCT__ 2124 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 2125 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2126 { 2127 PetscErrorCode ierr; 2128 PetscInt i,j,anzj; 2129 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2130 PetscInt an=A->cmap->N,am=A->rmap->N; 2131 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2132 2133 PetscFunctionBegin; 2134 /* Allocate space for symbolic transpose info and work array */ 2135 ierr = PetscCalloc1((an+1),&ati);CHKERRQ(ierr); 2136 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2137 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2138 2139 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2140 /* Note: offset by 1 for fast conversion into csr format. */ 2141 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2142 /* Form ati for csr format of A^T. */ 2143 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2144 2145 /* Copy ati into atfill so we have locations of the next free space in atj */ 2146 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2147 2148 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2149 for (i=0;i<am;i++) { 2150 anzj = ai[i+1] - ai[i]; 2151 for (j=0;j<anzj;j++) { 2152 atj[atfill[*aj]] = i; 2153 atfill[*aj++] += 1; 2154 } 2155 } 2156 2157 /* Clean up temporary space and complete requests. */ 2158 ierr = PetscFree(atfill);CHKERRQ(ierr); 2159 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2160 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2161 2162 b = (Mat_SeqAIJ*)((*B)->data); 2163 b->free_a = PETSC_FALSE; 2164 b->free_ij = PETSC_TRUE; 2165 b->nonew = 0; 2166 PetscFunctionReturn(0); 2167 } 2168 2169 #undef __FUNCT__ 2170 #define __FUNCT__ "MatTranspose_SeqAIJ" 2171 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2172 { 2173 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2174 Mat C; 2175 PetscErrorCode ierr; 2176 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2177 MatScalar *array = a->a; 2178 2179 PetscFunctionBegin; 2180 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 2181 2182 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2183 ierr = PetscCalloc1((1+A->cmap->n),&col);CHKERRQ(ierr); 2184 2185 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2186 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2187 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2188 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2189 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2190 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2191 ierr = PetscFree(col);CHKERRQ(ierr); 2192 } else { 2193 C = *B; 2194 } 2195 2196 for (i=0; i<m; i++) { 2197 len = ai[i+1]-ai[i]; 2198 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2199 array += len; 2200 aj += len; 2201 } 2202 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2203 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2204 2205 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2206 *B = C; 2207 } else { 2208 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 2209 } 2210 PetscFunctionReturn(0); 2211 } 2212 2213 #undef __FUNCT__ 2214 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2215 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2216 { 2217 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2218 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2219 MatScalar *va,*vb; 2220 PetscErrorCode ierr; 2221 PetscInt ma,na,mb,nb, i; 2222 2223 PetscFunctionBegin; 2224 bij = (Mat_SeqAIJ*) B->data; 2225 2226 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2227 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2228 if (ma!=nb || na!=mb) { 2229 *f = PETSC_FALSE; 2230 PetscFunctionReturn(0); 2231 } 2232 aii = aij->i; bii = bij->i; 2233 adx = aij->j; bdx = bij->j; 2234 va = aij->a; vb = bij->a; 2235 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2236 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2237 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2238 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2239 2240 *f = PETSC_TRUE; 2241 for (i=0; i<ma; i++) { 2242 while (aptr[i]<aii[i+1]) { 2243 PetscInt idc,idr; 2244 PetscScalar vc,vr; 2245 /* column/row index/value */ 2246 idc = adx[aptr[i]]; 2247 idr = bdx[bptr[idc]]; 2248 vc = va[aptr[i]]; 2249 vr = vb[bptr[idc]]; 2250 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2251 *f = PETSC_FALSE; 2252 goto done; 2253 } else { 2254 aptr[i]++; 2255 if (B || i!=idc) bptr[idc]++; 2256 } 2257 } 2258 } 2259 done: 2260 ierr = PetscFree(aptr);CHKERRQ(ierr); 2261 ierr = PetscFree(bptr);CHKERRQ(ierr); 2262 PetscFunctionReturn(0); 2263 } 2264 2265 #undef __FUNCT__ 2266 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2267 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2268 { 2269 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2270 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2271 MatScalar *va,*vb; 2272 PetscErrorCode ierr; 2273 PetscInt ma,na,mb,nb, i; 2274 2275 PetscFunctionBegin; 2276 bij = (Mat_SeqAIJ*) B->data; 2277 2278 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2279 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2280 if (ma!=nb || na!=mb) { 2281 *f = PETSC_FALSE; 2282 PetscFunctionReturn(0); 2283 } 2284 aii = aij->i; bii = bij->i; 2285 adx = aij->j; bdx = bij->j; 2286 va = aij->a; vb = bij->a; 2287 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2288 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2289 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2290 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2291 2292 *f = PETSC_TRUE; 2293 for (i=0; i<ma; i++) { 2294 while (aptr[i]<aii[i+1]) { 2295 PetscInt idc,idr; 2296 PetscScalar vc,vr; 2297 /* column/row index/value */ 2298 idc = adx[aptr[i]]; 2299 idr = bdx[bptr[idc]]; 2300 vc = va[aptr[i]]; 2301 vr = vb[bptr[idc]]; 2302 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2303 *f = PETSC_FALSE; 2304 goto done; 2305 } else { 2306 aptr[i]++; 2307 if (B || i!=idc) bptr[idc]++; 2308 } 2309 } 2310 } 2311 done: 2312 ierr = PetscFree(aptr);CHKERRQ(ierr); 2313 ierr = PetscFree(bptr);CHKERRQ(ierr); 2314 PetscFunctionReturn(0); 2315 } 2316 2317 #undef __FUNCT__ 2318 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2319 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2320 { 2321 PetscErrorCode ierr; 2322 2323 PetscFunctionBegin; 2324 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327 2328 #undef __FUNCT__ 2329 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2330 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2331 { 2332 PetscErrorCode ierr; 2333 2334 PetscFunctionBegin; 2335 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2336 PetscFunctionReturn(0); 2337 } 2338 2339 #undef __FUNCT__ 2340 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2341 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2342 { 2343 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2344 PetscScalar *l,*r,x; 2345 MatScalar *v; 2346 PetscErrorCode ierr; 2347 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2348 2349 PetscFunctionBegin; 2350 if (ll) { 2351 /* The local size is used so that VecMPI can be passed to this routine 2352 by MatDiagonalScale_MPIAIJ */ 2353 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2354 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2355 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2356 v = a->a; 2357 for (i=0; i<m; i++) { 2358 x = l[i]; 2359 M = a->i[i+1] - a->i[i]; 2360 for (j=0; j<M; j++) (*v++) *= x; 2361 } 2362 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2363 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2364 } 2365 if (rr) { 2366 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2367 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2368 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2369 v = a->a; jj = a->j; 2370 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2371 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2372 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2373 } 2374 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2375 PetscFunctionReturn(0); 2376 } 2377 2378 #undef __FUNCT__ 2379 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2380 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2381 { 2382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2383 PetscErrorCode ierr; 2384 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2385 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2386 const PetscInt *irow,*icol; 2387 PetscInt nrows,ncols; 2388 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2389 MatScalar *a_new,*mat_a; 2390 Mat C; 2391 PetscBool stride,sorted; 2392 2393 PetscFunctionBegin; 2394 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2395 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2396 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2397 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2398 2399 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2400 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2401 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2402 2403 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2404 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2405 if (stride && step == 1) { 2406 /* special case of contiguous rows */ 2407 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2408 /* loop over new rows determining lens and starting points */ 2409 for (i=0; i<nrows; i++) { 2410 kstart = ai[irow[i]]; 2411 kend = kstart + ailen[irow[i]]; 2412 for (k=kstart; k<kend; k++) { 2413 if (aj[k] >= first) { 2414 starts[i] = k; 2415 break; 2416 } 2417 } 2418 sum = 0; 2419 while (k < kend) { 2420 if (aj[k++] >= first+ncols) break; 2421 sum++; 2422 } 2423 lens[i] = sum; 2424 } 2425 /* create submatrix */ 2426 if (scall == MAT_REUSE_MATRIX) { 2427 PetscInt n_cols,n_rows; 2428 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2429 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2430 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2431 C = *B; 2432 } else { 2433 PetscInt rbs,cbs; 2434 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2435 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2436 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2437 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2438 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2439 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2440 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2441 } 2442 c = (Mat_SeqAIJ*)C->data; 2443 2444 /* loop over rows inserting into submatrix */ 2445 a_new = c->a; 2446 j_new = c->j; 2447 i_new = c->i; 2448 2449 for (i=0; i<nrows; i++) { 2450 ii = starts[i]; 2451 lensi = lens[i]; 2452 for (k=0; k<lensi; k++) { 2453 *j_new++ = aj[ii+k] - first; 2454 } 2455 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2456 a_new += lensi; 2457 i_new[i+1] = i_new[i] + lensi; 2458 c->ilen[i] = lensi; 2459 } 2460 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2461 } else { 2462 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2463 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2464 ierr = PetscMalloc1((1+nrows),&lens);CHKERRQ(ierr); 2465 for (i=0; i<ncols; i++) { 2466 #if defined(PETSC_USE_DEBUG) 2467 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); 2468 #endif 2469 smap[icol[i]] = i+1; 2470 } 2471 2472 /* determine lens of each row */ 2473 for (i=0; i<nrows; i++) { 2474 kstart = ai[irow[i]]; 2475 kend = kstart + a->ilen[irow[i]]; 2476 lens[i] = 0; 2477 for (k=kstart; k<kend; k++) { 2478 if (smap[aj[k]]) { 2479 lens[i]++; 2480 } 2481 } 2482 } 2483 /* Create and fill new matrix */ 2484 if (scall == MAT_REUSE_MATRIX) { 2485 PetscBool equal; 2486 2487 c = (Mat_SeqAIJ*)((*B)->data); 2488 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2489 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2490 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2491 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2492 C = *B; 2493 } else { 2494 PetscInt rbs,cbs; 2495 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2496 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2497 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2498 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2499 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2500 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2501 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2502 } 2503 c = (Mat_SeqAIJ*)(C->data); 2504 for (i=0; i<nrows; i++) { 2505 row = irow[i]; 2506 kstart = ai[row]; 2507 kend = kstart + a->ilen[row]; 2508 mat_i = c->i[i]; 2509 mat_j = c->j + mat_i; 2510 mat_a = c->a + mat_i; 2511 mat_ilen = c->ilen + i; 2512 for (k=kstart; k<kend; k++) { 2513 if ((tcol=smap[a->j[k]])) { 2514 *mat_j++ = tcol - 1; 2515 *mat_a++ = a->a[k]; 2516 (*mat_ilen)++; 2517 2518 } 2519 } 2520 } 2521 /* Free work space */ 2522 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2523 ierr = PetscFree(smap);CHKERRQ(ierr); 2524 ierr = PetscFree(lens);CHKERRQ(ierr); 2525 } 2526 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2527 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2528 2529 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2530 *B = C; 2531 PetscFunctionReturn(0); 2532 } 2533 2534 #undef __FUNCT__ 2535 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2536 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2537 { 2538 PetscErrorCode ierr; 2539 Mat B; 2540 2541 PetscFunctionBegin; 2542 if (scall == MAT_INITIAL_MATRIX) { 2543 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2544 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2545 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2546 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2547 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2548 *subMat = B; 2549 } else { 2550 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2551 } 2552 PetscFunctionReturn(0); 2553 } 2554 2555 #undef __FUNCT__ 2556 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2557 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2558 { 2559 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2560 PetscErrorCode ierr; 2561 Mat outA; 2562 PetscBool row_identity,col_identity; 2563 2564 PetscFunctionBegin; 2565 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2566 2567 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2568 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2569 2570 outA = inA; 2571 outA->factortype = MAT_FACTOR_LU; 2572 2573 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2574 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2575 2576 a->row = row; 2577 2578 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2579 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2580 2581 a->col = col; 2582 2583 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2584 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2585 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2586 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2587 2588 if (!a->solve_work) { /* this matrix may have been factored before */ 2589 ierr = PetscMalloc1((inA->rmap->n+1),&a->solve_work);CHKERRQ(ierr); 2590 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2591 } 2592 2593 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2594 if (row_identity && col_identity) { 2595 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2596 } else { 2597 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2598 } 2599 PetscFunctionReturn(0); 2600 } 2601 2602 #undef __FUNCT__ 2603 #define __FUNCT__ "MatScale_SeqAIJ" 2604 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2605 { 2606 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2607 PetscScalar oalpha = alpha; 2608 PetscErrorCode ierr; 2609 PetscBLASInt one = 1,bnz; 2610 2611 PetscFunctionBegin; 2612 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2613 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2614 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2615 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2616 PetscFunctionReturn(0); 2617 } 2618 2619 #undef __FUNCT__ 2620 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2621 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2622 { 2623 PetscErrorCode ierr; 2624 PetscInt i; 2625 2626 PetscFunctionBegin; 2627 if (scall == MAT_INITIAL_MATRIX) { 2628 ierr = PetscMalloc1((n+1),B);CHKERRQ(ierr); 2629 } 2630 2631 for (i=0; i<n; i++) { 2632 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2633 } 2634 PetscFunctionReturn(0); 2635 } 2636 2637 #undef __FUNCT__ 2638 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2639 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2640 { 2641 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2642 PetscErrorCode ierr; 2643 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2644 const PetscInt *idx; 2645 PetscInt start,end,*ai,*aj; 2646 PetscBT table; 2647 2648 PetscFunctionBegin; 2649 m = A->rmap->n; 2650 ai = a->i; 2651 aj = a->j; 2652 2653 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2654 2655 ierr = PetscMalloc1((m+1),&nidx);CHKERRQ(ierr); 2656 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2657 2658 for (i=0; i<is_max; i++) { 2659 /* Initialize the two local arrays */ 2660 isz = 0; 2661 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2662 2663 /* Extract the indices, assume there can be duplicate entries */ 2664 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2665 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2666 2667 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2668 for (j=0; j<n; ++j) { 2669 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2670 } 2671 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2672 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2673 2674 k = 0; 2675 for (j=0; j<ov; j++) { /* for each overlap */ 2676 n = isz; 2677 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2678 row = nidx[k]; 2679 start = ai[row]; 2680 end = ai[row+1]; 2681 for (l = start; l<end; l++) { 2682 val = aj[l]; 2683 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2684 } 2685 } 2686 } 2687 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2688 } 2689 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2690 ierr = PetscFree(nidx);CHKERRQ(ierr); 2691 PetscFunctionReturn(0); 2692 } 2693 2694 /* -------------------------------------------------------------- */ 2695 #undef __FUNCT__ 2696 #define __FUNCT__ "MatPermute_SeqAIJ" 2697 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2698 { 2699 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2700 PetscErrorCode ierr; 2701 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2702 const PetscInt *row,*col; 2703 PetscInt *cnew,j,*lens; 2704 IS icolp,irowp; 2705 PetscInt *cwork = NULL; 2706 PetscScalar *vwork = NULL; 2707 2708 PetscFunctionBegin; 2709 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2710 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2711 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2712 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2713 2714 /* determine lengths of permuted rows */ 2715 ierr = PetscMalloc1((m+1),&lens);CHKERRQ(ierr); 2716 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2717 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2718 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2719 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2720 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2721 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2722 ierr = PetscFree(lens);CHKERRQ(ierr); 2723 2724 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2725 for (i=0; i<m; i++) { 2726 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2727 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2728 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2729 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2730 } 2731 ierr = PetscFree(cnew);CHKERRQ(ierr); 2732 2733 (*B)->assembled = PETSC_FALSE; 2734 2735 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2736 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2737 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2738 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2739 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2740 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2741 PetscFunctionReturn(0); 2742 } 2743 2744 #undef __FUNCT__ 2745 #define __FUNCT__ "MatCopy_SeqAIJ" 2746 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2747 { 2748 PetscErrorCode ierr; 2749 2750 PetscFunctionBegin; 2751 /* If the two matrices have the same copy implementation, use fast copy. */ 2752 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2753 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2754 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2755 2756 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"); 2757 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2758 } else { 2759 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2760 } 2761 PetscFunctionReturn(0); 2762 } 2763 2764 #undef __FUNCT__ 2765 #define __FUNCT__ "MatSetUp_SeqAIJ" 2766 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2767 { 2768 PetscErrorCode ierr; 2769 2770 PetscFunctionBegin; 2771 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2772 PetscFunctionReturn(0); 2773 } 2774 2775 #undef __FUNCT__ 2776 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2777 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2778 { 2779 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2780 2781 PetscFunctionBegin; 2782 *array = a->a; 2783 PetscFunctionReturn(0); 2784 } 2785 2786 #undef __FUNCT__ 2787 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2788 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2789 { 2790 PetscFunctionBegin; 2791 PetscFunctionReturn(0); 2792 } 2793 2794 /* 2795 Computes the number of nonzeros per row needed for preallocation when X and Y 2796 have different nonzero structure. 2797 */ 2798 #undef __FUNCT__ 2799 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2800 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2801 { 2802 PetscInt i,m=Y->rmap->N; 2803 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2804 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2805 const PetscInt *xi = x->i,*yi = y->i; 2806 2807 PetscFunctionBegin; 2808 /* Set the number of nonzeros in the new matrix */ 2809 for (i=0; i<m; i++) { 2810 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2811 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2812 nnz[i] = 0; 2813 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2814 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2815 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2816 nnz[i]++; 2817 } 2818 for (; k<nzy; k++) nnz[i]++; 2819 } 2820 PetscFunctionReturn(0); 2821 } 2822 2823 #undef __FUNCT__ 2824 #define __FUNCT__ "MatAXPY_SeqAIJ" 2825 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2826 { 2827 PetscErrorCode ierr; 2828 PetscInt i; 2829 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2830 PetscBLASInt one=1,bnz; 2831 2832 PetscFunctionBegin; 2833 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2834 if (str == SAME_NONZERO_PATTERN) { 2835 PetscScalar alpha = a; 2836 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2837 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2838 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2839 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2840 if (y->xtoy && y->XtoY != X) { 2841 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2842 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2843 } 2844 if (!y->xtoy) { /* get xtoy */ 2845 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,NULL, y->i,y->j,NULL, &y->xtoy);CHKERRQ(ierr); 2846 y->XtoY = X; 2847 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2848 } 2849 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2850 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2851 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %D/%D = %g\n",x->nz,y->nz,(double)((PetscReal)(x->nz)/(y->nz+1)));CHKERRQ(ierr); 2852 } else { 2853 Mat B; 2854 PetscInt *nnz; 2855 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2856 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2857 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2858 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2859 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2860 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2861 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2862 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2863 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2864 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2865 ierr = PetscFree(nnz);CHKERRQ(ierr); 2866 } 2867 PetscFunctionReturn(0); 2868 } 2869 2870 #undef __FUNCT__ 2871 #define __FUNCT__ "MatConjugate_SeqAIJ" 2872 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2873 { 2874 #if defined(PETSC_USE_COMPLEX) 2875 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2876 PetscInt i,nz; 2877 PetscScalar *a; 2878 2879 PetscFunctionBegin; 2880 nz = aij->nz; 2881 a = aij->a; 2882 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2883 #else 2884 PetscFunctionBegin; 2885 #endif 2886 PetscFunctionReturn(0); 2887 } 2888 2889 #undef __FUNCT__ 2890 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2891 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2892 { 2893 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2894 PetscErrorCode ierr; 2895 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2896 PetscReal atmp; 2897 PetscScalar *x; 2898 MatScalar *aa; 2899 2900 PetscFunctionBegin; 2901 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2902 aa = a->a; 2903 ai = a->i; 2904 aj = a->j; 2905 2906 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2907 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2908 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2909 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2910 for (i=0; i<m; i++) { 2911 ncols = ai[1] - ai[0]; ai++; 2912 x[i] = 0.0; 2913 for (j=0; j<ncols; j++) { 2914 atmp = PetscAbsScalar(*aa); 2915 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2916 aa++; aj++; 2917 } 2918 } 2919 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2920 PetscFunctionReturn(0); 2921 } 2922 2923 #undef __FUNCT__ 2924 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2925 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2926 { 2927 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2928 PetscErrorCode ierr; 2929 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2930 PetscScalar *x; 2931 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 maximum is 0.0 or higher */ 2948 x[i] = 0.0; 2949 if (idx) { 2950 idx[i] = 0; /* in case ncols is zero */ 2951 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 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 #undef __FUNCT__ 2969 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2970 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2971 { 2972 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2973 PetscErrorCode ierr; 2974 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2975 PetscReal atmp; 2976 PetscScalar *x; 2977 MatScalar *aa; 2978 2979 PetscFunctionBegin; 2980 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2981 aa = a->a; 2982 ai = a->i; 2983 aj = a->j; 2984 2985 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2986 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2987 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2988 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); 2989 for (i=0; i<m; i++) { 2990 ncols = ai[1] - ai[0]; ai++; 2991 if (ncols) { 2992 /* Get first nonzero */ 2993 for (j = 0; j < ncols; j++) { 2994 atmp = PetscAbsScalar(aa[j]); 2995 if (atmp > 1.0e-12) { 2996 x[i] = atmp; 2997 if (idx) idx[i] = aj[j]; 2998 break; 2999 } 3000 } 3001 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3002 } else { 3003 x[i] = 0.0; if (idx) idx[i] = 0; 3004 } 3005 for (j = 0; j < ncols; j++) { 3006 atmp = PetscAbsScalar(*aa); 3007 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3008 aa++; aj++; 3009 } 3010 } 3011 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3012 PetscFunctionReturn(0); 3013 } 3014 3015 #undef __FUNCT__ 3016 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 3017 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3018 { 3019 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3020 PetscErrorCode ierr; 3021 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3022 PetscScalar *x; 3023 MatScalar *aa; 3024 3025 PetscFunctionBegin; 3026 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3027 aa = a->a; 3028 ai = a->i; 3029 aj = a->j; 3030 3031 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3032 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3033 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3034 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3035 for (i=0; i<m; i++) { 3036 ncols = ai[1] - ai[0]; ai++; 3037 if (ncols == A->cmap->n) { /* row is dense */ 3038 x[i] = *aa; if (idx) idx[i] = 0; 3039 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3040 x[i] = 0.0; 3041 if (idx) { /* find first implicit 0.0 in the row */ 3042 idx[i] = 0; /* in case ncols is zero */ 3043 for (j=0; j<ncols; j++) { 3044 if (aj[j] > j) { 3045 idx[i] = j; 3046 break; 3047 } 3048 } 3049 } 3050 } 3051 for (j=0; j<ncols; j++) { 3052 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3053 aa++; aj++; 3054 } 3055 } 3056 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3057 PetscFunctionReturn(0); 3058 } 3059 3060 #include <petscblaslapack.h> 3061 #include <petsc-private/kernels/blockinvert.h> 3062 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 3065 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3066 { 3067 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3068 PetscErrorCode ierr; 3069 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3070 MatScalar *diag,work[25],*v_work; 3071 PetscReal shift = 0.0; 3072 3073 PetscFunctionBegin; 3074 if (a->ibdiagvalid) { 3075 if (values) *values = a->ibdiag; 3076 PetscFunctionReturn(0); 3077 } 3078 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3079 if (!a->ibdiag) { 3080 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3081 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3082 } 3083 diag = a->ibdiag; 3084 if (values) *values = a->ibdiag; 3085 /* factor and invert each block */ 3086 switch (bs) { 3087 case 1: 3088 for (i=0; i<mbs; i++) { 3089 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3090 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3091 } 3092 break; 3093 case 2: 3094 for (i=0; i<mbs; i++) { 3095 ij[0] = 2*i; ij[1] = 2*i + 1; 3096 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3097 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 3098 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3099 diag += 4; 3100 } 3101 break; 3102 case 3: 3103 for (i=0; i<mbs; i++) { 3104 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3105 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3106 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 3107 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3108 diag += 9; 3109 } 3110 break; 3111 case 4: 3112 for (i=0; i<mbs; i++) { 3113 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3114 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3115 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 3116 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3117 diag += 16; 3118 } 3119 break; 3120 case 5: 3121 for (i=0; i<mbs; i++) { 3122 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3123 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3124 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 3125 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3126 diag += 25; 3127 } 3128 break; 3129 case 6: 3130 for (i=0; i<mbs; i++) { 3131 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; 3132 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3133 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 3134 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3135 diag += 36; 3136 } 3137 break; 3138 case 7: 3139 for (i=0; i<mbs; i++) { 3140 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3141 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3142 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3143 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3144 diag += 49; 3145 } 3146 break; 3147 default: 3148 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3149 for (i=0; i<mbs; i++) { 3150 for (j=0; j<bs; j++) { 3151 IJ[j] = bs*i + j; 3152 } 3153 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3154 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3155 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3156 diag += bs2; 3157 } 3158 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3159 } 3160 a->ibdiagvalid = PETSC_TRUE; 3161 PetscFunctionReturn(0); 3162 } 3163 3164 #undef __FUNCT__ 3165 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3166 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3167 { 3168 PetscErrorCode ierr; 3169 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3170 PetscScalar a; 3171 PetscInt m,n,i,j,col; 3172 3173 PetscFunctionBegin; 3174 if (!x->assembled) { 3175 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3176 for (i=0; i<m; i++) { 3177 for (j=0; j<aij->imax[i]; j++) { 3178 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3179 col = (PetscInt)(n*PetscRealPart(a)); 3180 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3181 } 3182 } 3183 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3184 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3185 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3186 PetscFunctionReturn(0); 3187 } 3188 3189 /* -------------------------------------------------------------------*/ 3190 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3191 MatGetRow_SeqAIJ, 3192 MatRestoreRow_SeqAIJ, 3193 MatMult_SeqAIJ, 3194 /* 4*/ MatMultAdd_SeqAIJ, 3195 MatMultTranspose_SeqAIJ, 3196 MatMultTransposeAdd_SeqAIJ, 3197 0, 3198 0, 3199 0, 3200 /* 10*/ 0, 3201 MatLUFactor_SeqAIJ, 3202 0, 3203 MatSOR_SeqAIJ, 3204 MatTranspose_SeqAIJ, 3205 /*1 5*/ MatGetInfo_SeqAIJ, 3206 MatEqual_SeqAIJ, 3207 MatGetDiagonal_SeqAIJ, 3208 MatDiagonalScale_SeqAIJ, 3209 MatNorm_SeqAIJ, 3210 /* 20*/ 0, 3211 MatAssemblyEnd_SeqAIJ, 3212 MatSetOption_SeqAIJ, 3213 MatZeroEntries_SeqAIJ, 3214 /* 24*/ MatZeroRows_SeqAIJ, 3215 0, 3216 0, 3217 0, 3218 0, 3219 /* 29*/ MatSetUp_SeqAIJ, 3220 0, 3221 0, 3222 0, 3223 0, 3224 /* 34*/ MatDuplicate_SeqAIJ, 3225 0, 3226 0, 3227 MatILUFactor_SeqAIJ, 3228 0, 3229 /* 39*/ MatAXPY_SeqAIJ, 3230 MatGetSubMatrices_SeqAIJ, 3231 MatIncreaseOverlap_SeqAIJ, 3232 MatGetValues_SeqAIJ, 3233 MatCopy_SeqAIJ, 3234 /* 44*/ MatGetRowMax_SeqAIJ, 3235 MatScale_SeqAIJ, 3236 0, 3237 MatDiagonalSet_SeqAIJ, 3238 MatZeroRowsColumns_SeqAIJ, 3239 /* 49*/ MatSetRandom_SeqAIJ, 3240 MatGetRowIJ_SeqAIJ, 3241 MatRestoreRowIJ_SeqAIJ, 3242 MatGetColumnIJ_SeqAIJ, 3243 MatRestoreColumnIJ_SeqAIJ, 3244 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3245 0, 3246 0, 3247 MatPermute_SeqAIJ, 3248 0, 3249 /* 59*/ 0, 3250 MatDestroy_SeqAIJ, 3251 MatView_SeqAIJ, 3252 0, 3253 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3254 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3255 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3256 0, 3257 0, 3258 0, 3259 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3260 MatGetRowMinAbs_SeqAIJ, 3261 0, 3262 MatSetColoring_SeqAIJ, 3263 0, 3264 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3265 MatFDColoringApply_AIJ, 3266 0, 3267 0, 3268 0, 3269 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3270 0, 3271 0, 3272 0, 3273 MatLoad_SeqAIJ, 3274 /* 84*/ MatIsSymmetric_SeqAIJ, 3275 MatIsHermitian_SeqAIJ, 3276 0, 3277 0, 3278 0, 3279 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3280 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3281 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3282 MatPtAP_SeqAIJ_SeqAIJ, 3283 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3284 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3285 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3286 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3287 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3288 0, 3289 /* 99*/ 0, 3290 0, 3291 0, 3292 MatConjugate_SeqAIJ, 3293 0, 3294 /*104*/ MatSetValuesRow_SeqAIJ, 3295 MatRealPart_SeqAIJ, 3296 MatImaginaryPart_SeqAIJ, 3297 0, 3298 0, 3299 /*109*/ MatMatSolve_SeqAIJ, 3300 0, 3301 MatGetRowMin_SeqAIJ, 3302 0, 3303 MatMissingDiagonal_SeqAIJ, 3304 /*114*/ 0, 3305 0, 3306 0, 3307 0, 3308 0, 3309 /*119*/ 0, 3310 0, 3311 0, 3312 0, 3313 MatGetMultiProcBlock_SeqAIJ, 3314 /*124*/ MatFindNonzeroRows_SeqAIJ, 3315 MatGetColumnNorms_SeqAIJ, 3316 MatInvertBlockDiagonal_SeqAIJ, 3317 0, 3318 0, 3319 /*129*/ 0, 3320 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3321 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3322 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3323 MatTransposeColoringCreate_SeqAIJ, 3324 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3325 MatTransColoringApplyDenToSp_SeqAIJ, 3326 MatRARt_SeqAIJ_SeqAIJ, 3327 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3328 MatRARtNumeric_SeqAIJ_SeqAIJ, 3329 /*139*/0, 3330 0, 3331 0, 3332 MatFDColoringSetUp_SeqXAIJ, 3333 MatFindOffBlockDiagonalEntries_SeqAIJ 3334 }; 3335 3336 #undef __FUNCT__ 3337 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3338 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3339 { 3340 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3341 PetscInt i,nz,n; 3342 3343 PetscFunctionBegin; 3344 nz = aij->maxnz; 3345 n = mat->rmap->n; 3346 for (i=0; i<nz; i++) { 3347 aij->j[i] = indices[i]; 3348 } 3349 aij->nz = nz; 3350 for (i=0; i<n; i++) { 3351 aij->ilen[i] = aij->imax[i]; 3352 } 3353 PetscFunctionReturn(0); 3354 } 3355 3356 #undef __FUNCT__ 3357 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3358 /*@ 3359 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3360 in the matrix. 3361 3362 Input Parameters: 3363 + mat - the SeqAIJ matrix 3364 - indices - the column indices 3365 3366 Level: advanced 3367 3368 Notes: 3369 This can be called if you have precomputed the nonzero structure of the 3370 matrix and want to provide it to the matrix object to improve the performance 3371 of the MatSetValues() operation. 3372 3373 You MUST have set the correct numbers of nonzeros per row in the call to 3374 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3375 3376 MUST be called before any calls to MatSetValues(); 3377 3378 The indices should start with zero, not one. 3379 3380 @*/ 3381 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3382 { 3383 PetscErrorCode ierr; 3384 3385 PetscFunctionBegin; 3386 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3387 PetscValidPointer(indices,2); 3388 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3389 PetscFunctionReturn(0); 3390 } 3391 3392 /* ----------------------------------------------------------------------------------------*/ 3393 3394 #undef __FUNCT__ 3395 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3396 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3397 { 3398 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3399 PetscErrorCode ierr; 3400 size_t nz = aij->i[mat->rmap->n]; 3401 3402 PetscFunctionBegin; 3403 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3404 3405 /* allocate space for values if not already there */ 3406 if (!aij->saved_values) { 3407 ierr = PetscMalloc1((nz+1),&aij->saved_values);CHKERRQ(ierr); 3408 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3409 } 3410 3411 /* copy values over */ 3412 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3413 PetscFunctionReturn(0); 3414 } 3415 3416 #undef __FUNCT__ 3417 #define __FUNCT__ "MatStoreValues" 3418 /*@ 3419 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3420 example, reuse of the linear part of a Jacobian, while recomputing the 3421 nonlinear portion. 3422 3423 Collect on Mat 3424 3425 Input Parameters: 3426 . mat - the matrix (currently only AIJ matrices support this option) 3427 3428 Level: advanced 3429 3430 Common Usage, with SNESSolve(): 3431 $ Create Jacobian matrix 3432 $ Set linear terms into matrix 3433 $ Apply boundary conditions to matrix, at this time matrix must have 3434 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3435 $ boundary conditions again will not change the nonzero structure 3436 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3437 $ ierr = MatStoreValues(mat); 3438 $ Call SNESSetJacobian() with matrix 3439 $ In your Jacobian routine 3440 $ ierr = MatRetrieveValues(mat); 3441 $ Set nonlinear terms in matrix 3442 3443 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3444 $ // build linear portion of Jacobian 3445 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3446 $ ierr = MatStoreValues(mat); 3447 $ loop over nonlinear iterations 3448 $ ierr = MatRetrieveValues(mat); 3449 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3450 $ // call MatAssemblyBegin/End() on matrix 3451 $ Solve linear system with Jacobian 3452 $ endloop 3453 3454 Notes: 3455 Matrix must already be assemblied before calling this routine 3456 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3457 calling this routine. 3458 3459 When this is called multiple times it overwrites the previous set of stored values 3460 and does not allocated additional space. 3461 3462 .seealso: MatRetrieveValues() 3463 3464 @*/ 3465 PetscErrorCode MatStoreValues(Mat mat) 3466 { 3467 PetscErrorCode ierr; 3468 3469 PetscFunctionBegin; 3470 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3471 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3472 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3473 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3474 PetscFunctionReturn(0); 3475 } 3476 3477 #undef __FUNCT__ 3478 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3479 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3480 { 3481 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3482 PetscErrorCode ierr; 3483 PetscInt nz = aij->i[mat->rmap->n]; 3484 3485 PetscFunctionBegin; 3486 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3487 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3488 /* copy values over */ 3489 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3490 PetscFunctionReturn(0); 3491 } 3492 3493 #undef __FUNCT__ 3494 #define __FUNCT__ "MatRetrieveValues" 3495 /*@ 3496 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3497 example, reuse of the linear part of a Jacobian, while recomputing the 3498 nonlinear portion. 3499 3500 Collect on Mat 3501 3502 Input Parameters: 3503 . mat - the matrix (currently on AIJ matrices support this option) 3504 3505 Level: advanced 3506 3507 .seealso: MatStoreValues() 3508 3509 @*/ 3510 PetscErrorCode MatRetrieveValues(Mat mat) 3511 { 3512 PetscErrorCode ierr; 3513 3514 PetscFunctionBegin; 3515 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3516 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3517 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3518 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3519 PetscFunctionReturn(0); 3520 } 3521 3522 3523 /* --------------------------------------------------------------------------------*/ 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "MatCreateSeqAIJ" 3526 /*@C 3527 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3528 (the default parallel PETSc format). For good matrix assembly performance 3529 the user should preallocate the matrix storage by setting the parameter nz 3530 (or the array nnz). By setting these parameters accurately, performance 3531 during matrix assembly can be increased by more than a factor of 50. 3532 3533 Collective on MPI_Comm 3534 3535 Input Parameters: 3536 + comm - MPI communicator, set to PETSC_COMM_SELF 3537 . m - number of rows 3538 . n - number of columns 3539 . nz - number of nonzeros per row (same for all rows) 3540 - nnz - array containing the number of nonzeros in the various rows 3541 (possibly different for each row) or NULL 3542 3543 Output Parameter: 3544 . A - the matrix 3545 3546 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3547 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3548 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3549 3550 Notes: 3551 If nnz is given then nz is ignored 3552 3553 The AIJ format (also called the Yale sparse matrix format or 3554 compressed row storage), is fully compatible with standard Fortran 77 3555 storage. That is, the stored row and column indices can begin at 3556 either one (as in Fortran) or zero. See the users' manual for details. 3557 3558 Specify the preallocated storage with either nz or nnz (not both). 3559 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3560 allocation. For large problems you MUST preallocate memory or you 3561 will get TERRIBLE performance, see the users' manual chapter on matrices. 3562 3563 By default, this format uses inodes (identical nodes) when possible, to 3564 improve numerical efficiency of matrix-vector products and solves. We 3565 search for consecutive rows with the same nonzero structure, thereby 3566 reusing matrix information to achieve increased efficiency. 3567 3568 Options Database Keys: 3569 + -mat_no_inode - Do not use inodes 3570 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3571 3572 Level: intermediate 3573 3574 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3575 3576 @*/ 3577 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3578 { 3579 PetscErrorCode ierr; 3580 3581 PetscFunctionBegin; 3582 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3583 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3584 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3585 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3586 PetscFunctionReturn(0); 3587 } 3588 3589 #undef __FUNCT__ 3590 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3591 /*@C 3592 MatSeqAIJSetPreallocation - For good matrix assembly performance 3593 the user should preallocate the matrix storage by setting the parameter nz 3594 (or the array nnz). By setting these parameters accurately, performance 3595 during matrix assembly can be increased by more than a factor of 50. 3596 3597 Collective on MPI_Comm 3598 3599 Input Parameters: 3600 + B - The matrix-free 3601 . nz - number of nonzeros per row (same for all rows) 3602 - nnz - array containing the number of nonzeros in the various rows 3603 (possibly different for each row) or NULL 3604 3605 Notes: 3606 If nnz is given then nz is ignored 3607 3608 The AIJ format (also called the Yale sparse matrix format or 3609 compressed row storage), is fully compatible with standard Fortran 77 3610 storage. That is, the stored row and column indices can begin at 3611 either one (as in Fortran) or zero. See the users' manual for details. 3612 3613 Specify the preallocated storage with either nz or nnz (not both). 3614 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3615 allocation. For large problems you MUST preallocate memory or you 3616 will get TERRIBLE performance, see the users' manual chapter on matrices. 3617 3618 You can call MatGetInfo() to get information on how effective the preallocation was; 3619 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3620 You can also run with the option -info and look for messages with the string 3621 malloc in them to see if additional memory allocation was needed. 3622 3623 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3624 entries or columns indices 3625 3626 By default, this format uses inodes (identical nodes) when possible, to 3627 improve numerical efficiency of matrix-vector products and solves. We 3628 search for consecutive rows with the same nonzero structure, thereby 3629 reusing matrix information to achieve increased efficiency. 3630 3631 Options Database Keys: 3632 + -mat_no_inode - Do not use inodes 3633 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3634 - -mat_aij_oneindex - Internally use indexing starting at 1 3635 rather than 0. Note that when calling MatSetValues(), 3636 the user still MUST index entries starting at 0! 3637 3638 Level: intermediate 3639 3640 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3641 3642 @*/ 3643 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3644 { 3645 PetscErrorCode ierr; 3646 3647 PetscFunctionBegin; 3648 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3649 PetscValidType(B,1); 3650 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3651 PetscFunctionReturn(0); 3652 } 3653 3654 #undef __FUNCT__ 3655 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3656 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3657 { 3658 Mat_SeqAIJ *b; 3659 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3660 PetscErrorCode ierr; 3661 PetscInt i; 3662 3663 PetscFunctionBegin; 3664 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3665 if (nz == MAT_SKIP_ALLOCATION) { 3666 skipallocation = PETSC_TRUE; 3667 nz = 0; 3668 } 3669 3670 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3671 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3672 3673 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3674 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3675 if (nnz) { 3676 for (i=0; i<B->rmap->n; i++) { 3677 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]); 3678 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); 3679 } 3680 } 3681 3682 B->preallocated = PETSC_TRUE; 3683 3684 b = (Mat_SeqAIJ*)B->data; 3685 3686 if (!skipallocation) { 3687 if (!b->imax) { 3688 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3689 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3690 } 3691 if (!nnz) { 3692 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3693 else if (nz < 0) nz = 1; 3694 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3695 nz = nz*B->rmap->n; 3696 } else { 3697 nz = 0; 3698 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3699 } 3700 /* b->ilen will count nonzeros in each row so far. */ 3701 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3702 3703 /* allocate the matrix space */ 3704 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3705 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3706 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3707 b->i[0] = 0; 3708 for (i=1; i<B->rmap->n+1; i++) { 3709 b->i[i] = b->i[i-1] + b->imax[i-1]; 3710 } 3711 b->singlemalloc = PETSC_TRUE; 3712 b->free_a = PETSC_TRUE; 3713 b->free_ij = PETSC_TRUE; 3714 #if defined(PETSC_THREADCOMM_ACTIVE) 3715 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3716 #endif 3717 } else { 3718 b->free_a = PETSC_FALSE; 3719 b->free_ij = PETSC_FALSE; 3720 } 3721 3722 b->nz = 0; 3723 b->maxnz = nz; 3724 B->info.nz_unneeded = (double)b->maxnz; 3725 if (realalloc) { 3726 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3727 } 3728 PetscFunctionReturn(0); 3729 } 3730 3731 #undef __FUNCT__ 3732 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3733 /*@ 3734 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3735 3736 Input Parameters: 3737 + B - the matrix 3738 . i - the indices into j for the start of each row (starts with zero) 3739 . j - the column indices for each row (starts with zero) these must be sorted for each row 3740 - v - optional values in the matrix 3741 3742 Level: developer 3743 3744 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3745 3746 .keywords: matrix, aij, compressed row, sparse, sequential 3747 3748 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3749 @*/ 3750 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3751 { 3752 PetscErrorCode ierr; 3753 3754 PetscFunctionBegin; 3755 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3756 PetscValidType(B,1); 3757 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3758 PetscFunctionReturn(0); 3759 } 3760 3761 #undef __FUNCT__ 3762 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3763 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3764 { 3765 PetscInt i; 3766 PetscInt m,n; 3767 PetscInt nz; 3768 PetscInt *nnz, nz_max = 0; 3769 PetscScalar *values; 3770 PetscErrorCode ierr; 3771 3772 PetscFunctionBegin; 3773 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3774 3775 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3776 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3777 3778 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3779 ierr = PetscMalloc1((m+1), &nnz);CHKERRQ(ierr); 3780 for (i = 0; i < m; i++) { 3781 nz = Ii[i+1]- Ii[i]; 3782 nz_max = PetscMax(nz_max, nz); 3783 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3784 nnz[i] = nz; 3785 } 3786 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3787 ierr = PetscFree(nnz);CHKERRQ(ierr); 3788 3789 if (v) { 3790 values = (PetscScalar*) v; 3791 } else { 3792 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3793 } 3794 3795 for (i = 0; i < m; i++) { 3796 nz = Ii[i+1] - Ii[i]; 3797 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3798 } 3799 3800 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3801 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3802 3803 if (!v) { 3804 ierr = PetscFree(values);CHKERRQ(ierr); 3805 } 3806 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3807 PetscFunctionReturn(0); 3808 } 3809 3810 #include <../src/mat/impls/dense/seq/dense.h> 3811 #include <petsc-private/kernels/petscaxpy.h> 3812 3813 #undef __FUNCT__ 3814 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3815 /* 3816 Computes (B'*A')' since computing B*A directly is untenable 3817 3818 n p p 3819 ( ) ( ) ( ) 3820 m ( A ) * n ( B ) = m ( C ) 3821 ( ) ( ) ( ) 3822 3823 */ 3824 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3825 { 3826 PetscErrorCode ierr; 3827 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3828 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3829 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3830 PetscInt i,n,m,q,p; 3831 const PetscInt *ii,*idx; 3832 const PetscScalar *b,*a,*a_q; 3833 PetscScalar *c,*c_q; 3834 3835 PetscFunctionBegin; 3836 m = A->rmap->n; 3837 n = A->cmap->n; 3838 p = B->cmap->n; 3839 a = sub_a->v; 3840 b = sub_b->a; 3841 c = sub_c->v; 3842 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3843 3844 ii = sub_b->i; 3845 idx = sub_b->j; 3846 for (i=0; i<n; i++) { 3847 q = ii[i+1] - ii[i]; 3848 while (q-->0) { 3849 c_q = c + m*(*idx); 3850 a_q = a + m*i; 3851 PetscKernelAXPY(c_q,*b,a_q,m); 3852 idx++; 3853 b++; 3854 } 3855 } 3856 PetscFunctionReturn(0); 3857 } 3858 3859 #undef __FUNCT__ 3860 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3861 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3862 { 3863 PetscErrorCode ierr; 3864 PetscInt m=A->rmap->n,n=B->cmap->n; 3865 Mat Cmat; 3866 3867 PetscFunctionBegin; 3868 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); 3869 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3870 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3871 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3872 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3873 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3874 3875 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3876 3877 *C = Cmat; 3878 PetscFunctionReturn(0); 3879 } 3880 3881 /* ----------------------------------------------------------------*/ 3882 #undef __FUNCT__ 3883 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3884 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3885 { 3886 PetscErrorCode ierr; 3887 3888 PetscFunctionBegin; 3889 if (scall == MAT_INITIAL_MATRIX) { 3890 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3891 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3892 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3893 } 3894 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3895 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3896 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3897 PetscFunctionReturn(0); 3898 } 3899 3900 3901 /*MC 3902 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3903 based on compressed sparse row format. 3904 3905 Options Database Keys: 3906 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3907 3908 Level: beginner 3909 3910 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3911 M*/ 3912 3913 /*MC 3914 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3915 3916 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3917 and MATMPIAIJ otherwise. As a result, for single process communicators, 3918 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3919 for communicators controlling multiple processes. It is recommended that you call both of 3920 the above preallocation routines for simplicity. 3921 3922 Options Database Keys: 3923 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3924 3925 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3926 enough exist. 3927 3928 Level: beginner 3929 3930 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3931 M*/ 3932 3933 /*MC 3934 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3935 3936 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3937 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3938 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3939 for communicators controlling multiple processes. It is recommended that you call both of 3940 the above preallocation routines for simplicity. 3941 3942 Options Database Keys: 3943 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3944 3945 Level: beginner 3946 3947 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3948 M*/ 3949 3950 #if defined(PETSC_HAVE_PASTIX) 3951 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3952 #endif 3953 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3954 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*); 3955 #endif 3956 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3957 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3958 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3959 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*); 3960 #if defined(PETSC_HAVE_MUMPS) 3961 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3962 #endif 3963 #if defined(PETSC_HAVE_SUPERLU) 3964 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3965 #endif 3966 #if defined(PETSC_HAVE_SUPERLU_DIST) 3967 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3968 #endif 3969 #if defined(PETSC_HAVE_SUITESPARSE) 3970 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3971 #endif 3972 #if defined(PETSC_HAVE_SUITESPARSE) 3973 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3974 #endif 3975 #if defined(PETSC_HAVE_SUITESPARSE) 3976 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*); 3977 #endif 3978 #if defined(PETSC_HAVE_LUSOL) 3979 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3980 #endif 3981 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3982 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3983 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3984 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3985 #endif 3986 #if defined(PETSC_HAVE_CLIQUE) 3987 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 3988 #endif 3989 3990 3991 #undef __FUNCT__ 3992 #define __FUNCT__ "MatSeqAIJGetArray" 3993 /*@C 3994 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3995 3996 Not Collective 3997 3998 Input Parameter: 3999 . mat - a MATSEQDENSE matrix 4000 4001 Output Parameter: 4002 . array - pointer to the data 4003 4004 Level: intermediate 4005 4006 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4007 @*/ 4008 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4009 { 4010 PetscErrorCode ierr; 4011 4012 PetscFunctionBegin; 4013 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4014 PetscFunctionReturn(0); 4015 } 4016 4017 #undef __FUNCT__ 4018 #define __FUNCT__ "MatSeqAIJRestoreArray" 4019 /*@C 4020 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 4021 4022 Not Collective 4023 4024 Input Parameters: 4025 . mat - a MATSEQDENSE matrix 4026 . array - pointer to the data 4027 4028 Level: intermediate 4029 4030 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4031 @*/ 4032 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4033 { 4034 PetscErrorCode ierr; 4035 4036 PetscFunctionBegin; 4037 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4038 PetscFunctionReturn(0); 4039 } 4040 4041 #undef __FUNCT__ 4042 #define __FUNCT__ "MatCreate_SeqAIJ" 4043 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4044 { 4045 Mat_SeqAIJ *b; 4046 PetscErrorCode ierr; 4047 PetscMPIInt size; 4048 4049 PetscFunctionBegin; 4050 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4051 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4052 4053 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4054 4055 B->data = (void*)b; 4056 4057 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4058 4059 b->row = 0; 4060 b->col = 0; 4061 b->icol = 0; 4062 b->reallocs = 0; 4063 b->ignorezeroentries = PETSC_FALSE; 4064 b->roworiented = PETSC_TRUE; 4065 b->nonew = 0; 4066 b->diag = 0; 4067 b->solve_work = 0; 4068 B->spptr = 0; 4069 b->saved_values = 0; 4070 b->idiag = 0; 4071 b->mdiag = 0; 4072 b->ssor_work = 0; 4073 b->omega = 1.0; 4074 b->fshift = 0.0; 4075 b->idiagvalid = PETSC_FALSE; 4076 b->ibdiagvalid = PETSC_FALSE; 4077 b->keepnonzeropattern = PETSC_FALSE; 4078 b->xtoy = 0; 4079 b->XtoY = 0; 4080 4081 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4082 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4083 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4084 4085 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4086 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 4087 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4088 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4089 #endif 4090 #if defined(PETSC_HAVE_PASTIX) 4091 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 4092 #endif 4093 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 4094 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 4095 #endif 4096 #if defined(PETSC_HAVE_SUPERLU) 4097 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 4098 #endif 4099 #if defined(PETSC_HAVE_SUPERLU_DIST) 4100 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 4101 #endif 4102 #if defined(PETSC_HAVE_MUMPS) 4103 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_aij_mumps);CHKERRQ(ierr); 4104 #endif 4105 #if defined(PETSC_HAVE_SUITESPARSE) 4106 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 4107 #endif 4108 #if defined(PETSC_HAVE_SUITESPARSE) 4109 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 4110 #endif 4111 #if defined(PETSC_HAVE_SUITESPARSE) 4112 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_klu_C",MatGetFactor_seqaij_klu);CHKERRQ(ierr); 4113 #endif 4114 #if defined(PETSC_HAVE_LUSOL) 4115 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 4116 #endif 4117 #if defined(PETSC_HAVE_CLIQUE) 4118 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C",MatGetFactor_aij_clique);CHKERRQ(ierr); 4119 #endif 4120 4121 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 4122 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 4123 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 4124 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4125 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4126 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4127 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4128 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4129 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4130 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4131 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4132 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4133 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4134 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4135 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4136 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4137 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4138 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4139 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4140 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4141 PetscFunctionReturn(0); 4142 } 4143 4144 #undef __FUNCT__ 4145 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4146 /* 4147 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4148 */ 4149 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4150 { 4151 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4152 PetscErrorCode ierr; 4153 PetscInt i,m = A->rmap->n; 4154 4155 PetscFunctionBegin; 4156 c = (Mat_SeqAIJ*)C->data; 4157 4158 C->factortype = A->factortype; 4159 c->row = 0; 4160 c->col = 0; 4161 c->icol = 0; 4162 c->reallocs = 0; 4163 4164 C->assembled = PETSC_TRUE; 4165 4166 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4167 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4168 4169 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4170 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4171 for (i=0; i<m; i++) { 4172 c->imax[i] = a->imax[i]; 4173 c->ilen[i] = a->ilen[i]; 4174 } 4175 4176 /* allocate the matrix space */ 4177 if (mallocmatspace) { 4178 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4179 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4180 4181 c->singlemalloc = PETSC_TRUE; 4182 4183 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4184 if (m > 0) { 4185 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4186 if (cpvalues == MAT_COPY_VALUES) { 4187 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4188 } else { 4189 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4190 } 4191 } 4192 } 4193 4194 c->ignorezeroentries = a->ignorezeroentries; 4195 c->roworiented = a->roworiented; 4196 c->nonew = a->nonew; 4197 if (a->diag) { 4198 ierr = PetscMalloc1((m+1),&c->diag);CHKERRQ(ierr); 4199 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4200 for (i=0; i<m; i++) { 4201 c->diag[i] = a->diag[i]; 4202 } 4203 } else c->diag = 0; 4204 4205 c->solve_work = 0; 4206 c->saved_values = 0; 4207 c->idiag = 0; 4208 c->ssor_work = 0; 4209 c->keepnonzeropattern = a->keepnonzeropattern; 4210 c->free_a = PETSC_TRUE; 4211 c->free_ij = PETSC_TRUE; 4212 c->xtoy = 0; 4213 c->XtoY = 0; 4214 4215 c->rmax = a->rmax; 4216 c->nz = a->nz; 4217 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4218 C->preallocated = PETSC_TRUE; 4219 4220 c->compressedrow.use = a->compressedrow.use; 4221 c->compressedrow.nrows = a->compressedrow.nrows; 4222 if (a->compressedrow.use) { 4223 i = a->compressedrow.nrows; 4224 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4225 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4226 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4227 } else { 4228 c->compressedrow.use = PETSC_FALSE; 4229 c->compressedrow.i = NULL; 4230 c->compressedrow.rindex = NULL; 4231 } 4232 C->nonzerostate = A->nonzerostate; 4233 4234 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4235 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4236 PetscFunctionReturn(0); 4237 } 4238 4239 #undef __FUNCT__ 4240 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4241 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4242 { 4243 PetscErrorCode ierr; 4244 4245 PetscFunctionBegin; 4246 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4247 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4248 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4249 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4250 } 4251 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4252 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4253 PetscFunctionReturn(0); 4254 } 4255 4256 #undef __FUNCT__ 4257 #define __FUNCT__ "MatLoad_SeqAIJ" 4258 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4259 { 4260 Mat_SeqAIJ *a; 4261 PetscErrorCode ierr; 4262 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4263 int fd; 4264 PetscMPIInt size; 4265 MPI_Comm comm; 4266 PetscInt bs = 1; 4267 4268 PetscFunctionBegin; 4269 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4270 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4271 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4272 4273 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4274 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4275 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4276 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4277 4278 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4279 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4280 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4281 M = header[1]; N = header[2]; nz = header[3]; 4282 4283 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4284 4285 /* read in row lengths */ 4286 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4287 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4288 4289 /* check if sum of rowlengths is same as nz */ 4290 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4291 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); 4292 4293 /* set global size if not set already*/ 4294 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4295 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4296 } else { 4297 /* if sizes and type are already set, check if the vector global sizes are correct */ 4298 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4299 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4300 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4301 } 4302 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); 4303 } 4304 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4305 a = (Mat_SeqAIJ*)newMat->data; 4306 4307 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4308 4309 /* read in nonzero values */ 4310 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4311 4312 /* set matrix "i" values */ 4313 a->i[0] = 0; 4314 for (i=1; i<= M; i++) { 4315 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4316 a->ilen[i-1] = rowlengths[i-1]; 4317 } 4318 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4319 4320 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4321 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4322 PetscFunctionReturn(0); 4323 } 4324 4325 #undef __FUNCT__ 4326 #define __FUNCT__ "MatEqual_SeqAIJ" 4327 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4328 { 4329 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4330 PetscErrorCode ierr; 4331 #if defined(PETSC_USE_COMPLEX) 4332 PetscInt k; 4333 #endif 4334 4335 PetscFunctionBegin; 4336 /* If the matrix dimensions are not equal,or no of nonzeros */ 4337 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4338 *flg = PETSC_FALSE; 4339 PetscFunctionReturn(0); 4340 } 4341 4342 /* if the a->i are the same */ 4343 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4344 if (!*flg) PetscFunctionReturn(0); 4345 4346 /* if a->j are the same */ 4347 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4348 if (!*flg) PetscFunctionReturn(0); 4349 4350 /* if a->a are the same */ 4351 #if defined(PETSC_USE_COMPLEX) 4352 for (k=0; k<a->nz; k++) { 4353 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4354 *flg = PETSC_FALSE; 4355 PetscFunctionReturn(0); 4356 } 4357 } 4358 #else 4359 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4360 #endif 4361 PetscFunctionReturn(0); 4362 } 4363 4364 #undef __FUNCT__ 4365 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4366 /*@ 4367 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4368 provided by the user. 4369 4370 Collective on MPI_Comm 4371 4372 Input Parameters: 4373 + comm - must be an MPI communicator of size 1 4374 . m - number of rows 4375 . n - number of columns 4376 . i - row indices 4377 . j - column indices 4378 - a - matrix values 4379 4380 Output Parameter: 4381 . mat - the matrix 4382 4383 Level: intermediate 4384 4385 Notes: 4386 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4387 once the matrix is destroyed and not before 4388 4389 You cannot set new nonzero locations into this matrix, that will generate an error. 4390 4391 The i and j indices are 0 based 4392 4393 The format which is used for the sparse matrix input, is equivalent to a 4394 row-major ordering.. i.e for the following matrix, the input data expected is 4395 as shown: 4396 4397 1 0 0 4398 2 0 3 4399 4 5 6 4400 4401 i = {0,1,3,6} [size = nrow+1 = 3+1] 4402 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4403 v = {1,2,3,4,5,6} [size = nz = 6] 4404 4405 4406 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4407 4408 @*/ 4409 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4410 { 4411 PetscErrorCode ierr; 4412 PetscInt ii; 4413 Mat_SeqAIJ *aij; 4414 #if defined(PETSC_USE_DEBUG) 4415 PetscInt jj; 4416 #endif 4417 4418 PetscFunctionBegin; 4419 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4420 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4421 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4422 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4423 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4424 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4425 aij = (Mat_SeqAIJ*)(*mat)->data; 4426 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4427 4428 aij->i = i; 4429 aij->j = j; 4430 aij->a = a; 4431 aij->singlemalloc = PETSC_FALSE; 4432 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4433 aij->free_a = PETSC_FALSE; 4434 aij->free_ij = PETSC_FALSE; 4435 4436 for (ii=0; ii<m; ii++) { 4437 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4438 #if defined(PETSC_USE_DEBUG) 4439 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]); 4440 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4441 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); 4442 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); 4443 } 4444 #endif 4445 } 4446 #if defined(PETSC_USE_DEBUG) 4447 for (ii=0; ii<aij->i[m]; ii++) { 4448 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4449 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]); 4450 } 4451 #endif 4452 4453 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4454 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4455 PetscFunctionReturn(0); 4456 } 4457 #undef __FUNCT__ 4458 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4459 /*@C 4460 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4461 provided by the user. 4462 4463 Collective on MPI_Comm 4464 4465 Input Parameters: 4466 + comm - must be an MPI communicator of size 1 4467 . m - number of rows 4468 . n - number of columns 4469 . i - row indices 4470 . j - column indices 4471 . a - matrix values 4472 . nz - number of nonzeros 4473 - idx - 0 or 1 based 4474 4475 Output Parameter: 4476 . mat - the matrix 4477 4478 Level: intermediate 4479 4480 Notes: 4481 The i and j indices are 0 based 4482 4483 The format which is used for the sparse matrix input, is equivalent to a 4484 row-major ordering.. i.e for the following matrix, the input data expected is 4485 as shown: 4486 4487 1 0 0 4488 2 0 3 4489 4 5 6 4490 4491 i = {0,1,1,2,2,2} 4492 j = {0,0,2,0,1,2} 4493 v = {1,2,3,4,5,6} 4494 4495 4496 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4497 4498 @*/ 4499 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4500 { 4501 PetscErrorCode ierr; 4502 PetscInt ii, *nnz, one = 1,row,col; 4503 4504 4505 PetscFunctionBegin; 4506 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4507 for (ii = 0; ii < nz; ii++) { 4508 nnz[i[ii] - !!idx] += 1; 4509 } 4510 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4511 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4512 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4513 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4514 for (ii = 0; ii < nz; ii++) { 4515 if (idx) { 4516 row = i[ii] - 1; 4517 col = j[ii] - 1; 4518 } else { 4519 row = i[ii]; 4520 col = j[ii]; 4521 } 4522 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4523 } 4524 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4525 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4526 ierr = PetscFree(nnz);CHKERRQ(ierr); 4527 PetscFunctionReturn(0); 4528 } 4529 4530 #undef __FUNCT__ 4531 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4532 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4533 { 4534 PetscErrorCode ierr; 4535 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4536 4537 PetscFunctionBegin; 4538 if (coloring->ctype == IS_COLORING_GLOBAL) { 4539 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4540 a->coloring = coloring; 4541 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4542 PetscInt i,*larray; 4543 ISColoring ocoloring; 4544 ISColoringValue *colors; 4545 4546 /* set coloring for diagonal portion */ 4547 ierr = PetscMalloc1(A->cmap->n,&larray);CHKERRQ(ierr); 4548 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4549 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4550 ierr = PetscMalloc1(A->cmap->n,&colors);CHKERRQ(ierr); 4551 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4552 ierr = PetscFree(larray);CHKERRQ(ierr); 4553 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4554 a->coloring = ocoloring; 4555 } 4556 PetscFunctionReturn(0); 4557 } 4558 4559 #undef __FUNCT__ 4560 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4561 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4562 { 4563 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4564 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4565 MatScalar *v = a->a; 4566 PetscScalar *values = (PetscScalar*)advalues; 4567 ISColoringValue *color; 4568 4569 PetscFunctionBegin; 4570 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4571 color = a->coloring->colors; 4572 /* loop over rows */ 4573 for (i=0; i<m; i++) { 4574 nz = ii[i+1] - ii[i]; 4575 /* loop over columns putting computed value into matrix */ 4576 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4577 values += nl; /* jump to next row of derivatives */ 4578 } 4579 PetscFunctionReturn(0); 4580 } 4581 4582 #undef __FUNCT__ 4583 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4584 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4585 { 4586 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4587 PetscErrorCode ierr; 4588 4589 PetscFunctionBegin; 4590 a->idiagvalid = PETSC_FALSE; 4591 a->ibdiagvalid = PETSC_FALSE; 4592 4593 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4594 PetscFunctionReturn(0); 4595 } 4596 4597 /* 4598 Special version for direct calls from Fortran 4599 */ 4600 #include <petsc-private/fortranimpl.h> 4601 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4602 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4603 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4604 #define matsetvaluesseqaij_ matsetvaluesseqaij 4605 #endif 4606 4607 /* Change these macros so can be used in void function */ 4608 #undef CHKERRQ 4609 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4610 #undef SETERRQ2 4611 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4612 #undef SETERRQ3 4613 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4614 4615 #undef __FUNCT__ 4616 #define __FUNCT__ "matsetvaluesseqaij_" 4617 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) 4618 { 4619 Mat A = *AA; 4620 PetscInt m = *mm, n = *nn; 4621 InsertMode is = *isis; 4622 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4623 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4624 PetscInt *imax,*ai,*ailen; 4625 PetscErrorCode ierr; 4626 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4627 MatScalar *ap,value,*aa; 4628 PetscBool ignorezeroentries = a->ignorezeroentries; 4629 PetscBool roworiented = a->roworiented; 4630 4631 PetscFunctionBegin; 4632 MatCheckPreallocated(A,1); 4633 imax = a->imax; 4634 ai = a->i; 4635 ailen = a->ilen; 4636 aj = a->j; 4637 aa = a->a; 4638 4639 for (k=0; k<m; k++) { /* loop over added rows */ 4640 row = im[k]; 4641 if (row < 0) continue; 4642 #if defined(PETSC_USE_DEBUG) 4643 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4644 #endif 4645 rp = aj + ai[row]; ap = aa + ai[row]; 4646 rmax = imax[row]; nrow = ailen[row]; 4647 low = 0; 4648 high = nrow; 4649 for (l=0; l<n; l++) { /* loop over added columns */ 4650 if (in[l] < 0) continue; 4651 #if defined(PETSC_USE_DEBUG) 4652 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4653 #endif 4654 col = in[l]; 4655 if (roworiented) value = v[l + k*n]; 4656 else value = v[k + l*m]; 4657 4658 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4659 4660 if (col <= lastcol) low = 0; 4661 else high = nrow; 4662 lastcol = col; 4663 while (high-low > 5) { 4664 t = (low+high)/2; 4665 if (rp[t] > col) high = t; 4666 else low = t; 4667 } 4668 for (i=low; i<high; i++) { 4669 if (rp[i] > col) break; 4670 if (rp[i] == col) { 4671 if (is == ADD_VALUES) ap[i] += value; 4672 else ap[i] = value; 4673 goto noinsert; 4674 } 4675 } 4676 if (value == 0.0 && ignorezeroentries) goto noinsert; 4677 if (nonew == 1) goto noinsert; 4678 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4679 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4680 N = nrow++ - 1; a->nz++; high++; 4681 /* shift up all the later entries in this row */ 4682 for (ii=N; ii>=i; ii--) { 4683 rp[ii+1] = rp[ii]; 4684 ap[ii+1] = ap[ii]; 4685 } 4686 rp[i] = col; 4687 ap[i] = value; 4688 A->nonzerostate++; 4689 noinsert:; 4690 low = i + 1; 4691 } 4692 ailen[row] = nrow; 4693 } 4694 PetscFunctionReturnVoid(); 4695 } 4696 4697 4698