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