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