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