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