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