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