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