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