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