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,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 212 } else if (oshift == 1) { 213 PetscInt *tia; 214 PetscInt nz = a->i[A->rmap->n]; 215 /* malloc space and add 1 to i and j indices */ 216 ierr = PetscMalloc1(A->rmap->n+1,&tia);CHKERRQ(ierr); 217 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 218 *ia = tia; 219 if (ja) { 220 PetscInt *tja; 221 ierr = PetscMalloc1(nz+1,&tja);CHKERRQ(ierr); 222 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 223 *ja = tja; 224 } 225 } else { 226 *ia = a->i; 227 if (ja) *ja = a->j; 228 } 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 234 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 if (!ia) PetscFunctionReturn(0); 240 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 241 ierr = PetscFree(*ia);CHKERRQ(ierr); 242 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 243 } 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNCT__ 248 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 249 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 250 { 251 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 252 PetscErrorCode ierr; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 261 } else { 262 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 263 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 264 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 ierr = PetscFree(collengths);CHKERRQ(ierr); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 291 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 292 { 293 PetscErrorCode ierr; 294 295 PetscFunctionBegin; 296 if (!ia) PetscFunctionReturn(0); 297 298 ierr = PetscFree(*ia);CHKERRQ(ierr); 299 ierr = PetscFree(*ja);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 /* 304 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 305 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 306 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 307 */ 308 #undef __FUNCT__ 309 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 310 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 311 { 312 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 313 PetscErrorCode ierr; 314 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 315 PetscInt nz = a->i[m],row,*jj,mr,col; 316 PetscInt *cspidx; 317 318 PetscFunctionBegin; 319 *nn = n; 320 if (!ia) PetscFunctionReturn(0); 321 322 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 323 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 324 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 325 ierr = PetscMalloc1(nz+1,&cspidx);CHKERRQ(ierr); 326 jj = a->j; 327 for (i=0; i<nz; i++) { 328 collengths[jj[i]]++; 329 } 330 cia[0] = oshift; 331 for (i=0; i<n; i++) { 332 cia[i+1] = cia[i] + collengths[i]; 333 } 334 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 335 jj = a->j; 336 for (row=0; row<m; row++) { 337 mr = a->i[row+1] - a->i[row]; 338 for (i=0; i<mr; i++) { 339 col = *jj++; 340 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 341 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 342 } 343 } 344 ierr = PetscFree(collengths);CHKERRQ(ierr); 345 *ia = cia; *ja = cja; 346 *spidx = cspidx; 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 352 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 353 { 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 ierr = MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);CHKERRQ(ierr); 358 ierr = PetscFree(*spidx);CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 364 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 365 { 366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 367 PetscInt *ai = a->i; 368 PetscErrorCode ierr; 369 370 PetscFunctionBegin; 371 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 /* 376 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 377 378 - a single row of values is set with each call 379 - no row or column indices are negative or (in error) larger than the number of rows or columns 380 - the values are always added to the matrix, not set 381 - no new locations are introduced in the nonzero structure of the matrix 382 383 This does NOT assume the global column indices are sorted 384 385 */ 386 387 #include <petsc/private/isimpl.h> 388 #undef __FUNCT__ 389 #define __FUNCT__ "MatSeqAIJSetValuesLocalFast" 390 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 391 { 392 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 393 PetscInt low,high,t,row,nrow,i,col,l; 394 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 395 PetscInt lastcol = -1; 396 MatScalar *ap,value,*aa = a->a; 397 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 398 399 row = ridx[im[0]]; 400 rp = aj + ai[row]; 401 ap = aa + ai[row]; 402 nrow = ailen[row]; 403 low = 0; 404 high = nrow; 405 for (l=0; l<n; l++) { /* loop over added columns */ 406 col = cidx[in[l]]; 407 value = v[l]; 408 409 if (col <= lastcol) low = 0; 410 else high = nrow; 411 lastcol = col; 412 while (high-low > 5) { 413 t = (low+high)/2; 414 if (rp[t] > col) high = t; 415 else low = t; 416 } 417 for (i=low; i<high; i++) { 418 if (rp[i] == col) { 419 ap[i] += value; 420 low = i + 1; 421 break; 422 } 423 } 424 } 425 return 0; 426 } 427 428 #undef __FUNCT__ 429 #define __FUNCT__ "MatSetValues_SeqAIJ" 430 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 431 { 432 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 433 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 434 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 435 PetscErrorCode ierr; 436 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 437 MatScalar *ap,value,*aa = a->a; 438 PetscBool ignorezeroentries = a->ignorezeroentries; 439 PetscBool roworiented = a->roworiented; 440 441 PetscFunctionBegin; 442 for (k=0; k<m; k++) { /* loop over added rows */ 443 row = im[k]; 444 if (row < 0) continue; 445 #if defined(PETSC_USE_DEBUG) 446 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 447 #endif 448 rp = aj + ai[row]; ap = aa + ai[row]; 449 rmax = imax[row]; nrow = ailen[row]; 450 low = 0; 451 high = nrow; 452 for (l=0; l<n; l++) { /* loop over added columns */ 453 if (in[l] < 0) continue; 454 #if defined(PETSC_USE_DEBUG) 455 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 456 #endif 457 col = in[l]; 458 if (roworiented) { 459 value = v[l + k*n]; 460 } else { 461 value = v[k + l*m]; 462 } 463 if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES)) continue; 464 465 if (col <= lastcol) low = 0; 466 else high = nrow; 467 lastcol = col; 468 while (high-low > 5) { 469 t = (low+high)/2; 470 if (rp[t] > col) high = t; 471 else low = t; 472 } 473 for (i=low; i<high; i++) { 474 if (rp[i] > col) break; 475 if (rp[i] == col) { 476 if (is == ADD_VALUES) ap[i] += value; 477 else ap[i] = value; 478 low = i + 1; 479 goto noinsert; 480 } 481 } 482 if (value == 0.0 && ignorezeroentries) goto noinsert; 483 if (nonew == 1) goto noinsert; 484 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 485 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 486 N = nrow++ - 1; a->nz++; high++; 487 /* shift up all the later entries in this row */ 488 for (ii=N; ii>=i; ii--) { 489 rp[ii+1] = rp[ii]; 490 ap[ii+1] = ap[ii]; 491 } 492 rp[i] = col; 493 ap[i] = value; 494 low = i + 1; 495 A->nonzerostate++; 496 noinsert:; 497 } 498 ailen[row] = nrow; 499 } 500 PetscFunctionReturn(0); 501 } 502 503 504 #undef __FUNCT__ 505 #define __FUNCT__ "MatGetValues_SeqAIJ" 506 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 507 { 508 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 509 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 510 PetscInt *ai = a->i,*ailen = a->ilen; 511 MatScalar *ap,*aa = a->a; 512 513 PetscFunctionBegin; 514 for (k=0; k<m; k++) { /* loop over rows */ 515 row = im[k]; 516 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 517 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 518 rp = aj + ai[row]; ap = aa + ai[row]; 519 nrow = ailen[row]; 520 for (l=0; l<n; l++) { /* loop over columns */ 521 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 522 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 523 col = in[l]; 524 high = nrow; low = 0; /* assume unsorted */ 525 while (high-low > 5) { 526 t = (low+high)/2; 527 if (rp[t] > col) high = t; 528 else low = t; 529 } 530 for (i=low; i<high; i++) { 531 if (rp[i] > col) break; 532 if (rp[i] == col) { 533 *v++ = ap[i]; 534 goto finished; 535 } 536 } 537 *v++ = 0.0; 538 finished:; 539 } 540 } 541 PetscFunctionReturn(0); 542 } 543 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatView_SeqAIJ_Binary" 547 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 548 { 549 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 550 PetscErrorCode ierr; 551 PetscInt i,*col_lens; 552 int fd; 553 FILE *file; 554 555 PetscFunctionBegin; 556 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 557 ierr = PetscMalloc1(4+A->rmap->n,&col_lens);CHKERRQ(ierr); 558 559 col_lens[0] = MAT_FILE_CLASSID; 560 col_lens[1] = A->rmap->n; 561 col_lens[2] = A->cmap->n; 562 col_lens[3] = a->nz; 563 564 /* store lengths of each row and write (including header) to file */ 565 for (i=0; i<A->rmap->n; i++) { 566 col_lens[4+i] = a->i[i+1] - a->i[i]; 567 } 568 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 569 ierr = PetscFree(col_lens);CHKERRQ(ierr); 570 571 /* store column indices (zero start index) */ 572 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 573 574 /* store nonzero values */ 575 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 576 577 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 578 if (file) { 579 fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs)); 580 } 581 PetscFunctionReturn(0); 582 } 583 584 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 588 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 589 { 590 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 591 PetscErrorCode ierr; 592 PetscInt i,j,m = A->rmap->n; 593 const char *name; 594 PetscViewerFormat format; 595 596 PetscFunctionBegin; 597 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 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1963 } else if (type == NORM_1) { 1964 PetscReal *tmp; 1965 PetscInt *jj = a->j; 1966 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 1967 *nrm = 0.0; 1968 for (j=0; j<a->nz; j++) { 1969 tmp[*jj++] += PetscAbsScalar(*v); v++; 1970 } 1971 for (j=0; j<A->cmap->n; j++) { 1972 if (tmp[j] > *nrm) *nrm = tmp[j]; 1973 } 1974 ierr = PetscFree(tmp);CHKERRQ(ierr); 1975 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1976 } else if (type == NORM_INFINITY) { 1977 *nrm = 0.0; 1978 for (j=0; j<A->rmap->n; j++) { 1979 v = a->a + a->i[j]; 1980 sum = 0.0; 1981 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1982 sum += PetscAbsScalar(*v); v++; 1983 } 1984 if (sum > *nrm) *nrm = sum; 1985 } 1986 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1987 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1988 PetscFunctionReturn(0); 1989 } 1990 1991 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 1994 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1995 { 1996 PetscErrorCode ierr; 1997 PetscInt i,j,anzj; 1998 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1999 PetscInt an=A->cmap->N,am=A->rmap->N; 2000 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2001 2002 PetscFunctionBegin; 2003 /* Allocate space for symbolic transpose info and work array */ 2004 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 2005 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2006 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2007 2008 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2009 /* Note: offset by 1 for fast conversion into csr format. */ 2010 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2011 /* Form ati for csr format of A^T. */ 2012 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2013 2014 /* Copy ati into atfill so we have locations of the next free space in atj */ 2015 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2016 2017 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2018 for (i=0;i<am;i++) { 2019 anzj = ai[i+1] - ai[i]; 2020 for (j=0;j<anzj;j++) { 2021 atj[atfill[*aj]] = i; 2022 atfill[*aj++] += 1; 2023 } 2024 } 2025 2026 /* Clean up temporary space and complete requests. */ 2027 ierr = PetscFree(atfill);CHKERRQ(ierr); 2028 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2029 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2030 2031 b = (Mat_SeqAIJ*)((*B)->data); 2032 b->free_a = PETSC_FALSE; 2033 b->free_ij = PETSC_TRUE; 2034 b->nonew = 0; 2035 PetscFunctionReturn(0); 2036 } 2037 2038 #undef __FUNCT__ 2039 #define __FUNCT__ "MatTranspose_SeqAIJ" 2040 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2041 { 2042 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2043 Mat C; 2044 PetscErrorCode ierr; 2045 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2046 MatScalar *array = a->a; 2047 2048 PetscFunctionBegin; 2049 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"); 2050 2051 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2052 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2053 2054 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2055 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2056 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2057 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2058 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2059 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2060 ierr = PetscFree(col);CHKERRQ(ierr); 2061 } else { 2062 C = *B; 2063 } 2064 2065 for (i=0; i<m; i++) { 2066 len = ai[i+1]-ai[i]; 2067 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2068 array += len; 2069 aj += len; 2070 } 2071 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2072 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2073 2074 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2075 *B = C; 2076 } else { 2077 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2078 } 2079 PetscFunctionReturn(0); 2080 } 2081 2082 #undef __FUNCT__ 2083 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2084 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2085 { 2086 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2087 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2088 MatScalar *va,*vb; 2089 PetscErrorCode ierr; 2090 PetscInt ma,na,mb,nb, i; 2091 2092 PetscFunctionBegin; 2093 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2094 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2095 if (ma!=nb || na!=mb) { 2096 *f = PETSC_FALSE; 2097 PetscFunctionReturn(0); 2098 } 2099 aii = aij->i; bii = bij->i; 2100 adx = aij->j; bdx = bij->j; 2101 va = aij->a; vb = bij->a; 2102 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2103 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2104 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2105 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2106 2107 *f = PETSC_TRUE; 2108 for (i=0; i<ma; i++) { 2109 while (aptr[i]<aii[i+1]) { 2110 PetscInt idc,idr; 2111 PetscScalar vc,vr; 2112 /* column/row index/value */ 2113 idc = adx[aptr[i]]; 2114 idr = bdx[bptr[idc]]; 2115 vc = va[aptr[i]]; 2116 vr = vb[bptr[idc]]; 2117 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2118 *f = PETSC_FALSE; 2119 goto done; 2120 } else { 2121 aptr[i]++; 2122 if (B || i!=idc) bptr[idc]++; 2123 } 2124 } 2125 } 2126 done: 2127 ierr = PetscFree(aptr);CHKERRQ(ierr); 2128 ierr = PetscFree(bptr);CHKERRQ(ierr); 2129 PetscFunctionReturn(0); 2130 } 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2134 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2135 { 2136 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2137 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2138 MatScalar *va,*vb; 2139 PetscErrorCode ierr; 2140 PetscInt ma,na,mb,nb, i; 2141 2142 PetscFunctionBegin; 2143 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2144 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2145 if (ma!=nb || na!=mb) { 2146 *f = PETSC_FALSE; 2147 PetscFunctionReturn(0); 2148 } 2149 aii = aij->i; bii = bij->i; 2150 adx = aij->j; bdx = bij->j; 2151 va = aij->a; vb = bij->a; 2152 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2153 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2154 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2155 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2156 2157 *f = PETSC_TRUE; 2158 for (i=0; i<ma; i++) { 2159 while (aptr[i]<aii[i+1]) { 2160 PetscInt idc,idr; 2161 PetscScalar vc,vr; 2162 /* column/row index/value */ 2163 idc = adx[aptr[i]]; 2164 idr = bdx[bptr[idc]]; 2165 vc = va[aptr[i]]; 2166 vr = vb[bptr[idc]]; 2167 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2168 *f = PETSC_FALSE; 2169 goto done; 2170 } else { 2171 aptr[i]++; 2172 if (B || i!=idc) bptr[idc]++; 2173 } 2174 } 2175 } 2176 done: 2177 ierr = PetscFree(aptr);CHKERRQ(ierr); 2178 ierr = PetscFree(bptr);CHKERRQ(ierr); 2179 PetscFunctionReturn(0); 2180 } 2181 2182 #undef __FUNCT__ 2183 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2184 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2185 { 2186 PetscErrorCode ierr; 2187 2188 PetscFunctionBegin; 2189 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2190 PetscFunctionReturn(0); 2191 } 2192 2193 #undef __FUNCT__ 2194 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2195 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2196 { 2197 PetscErrorCode ierr; 2198 2199 PetscFunctionBegin; 2200 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2201 PetscFunctionReturn(0); 2202 } 2203 2204 #undef __FUNCT__ 2205 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2206 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2207 { 2208 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2209 PetscScalar *l,*r,x; 2210 MatScalar *v; 2211 PetscErrorCode ierr; 2212 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2213 2214 PetscFunctionBegin; 2215 if (ll) { 2216 /* The local size is used so that VecMPI can be passed to this routine 2217 by MatDiagonalScale_MPIAIJ */ 2218 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2219 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2220 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2221 v = a->a; 2222 for (i=0; i<m; i++) { 2223 x = l[i]; 2224 M = a->i[i+1] - a->i[i]; 2225 for (j=0; j<M; j++) (*v++) *= x; 2226 } 2227 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2228 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2229 } 2230 if (rr) { 2231 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2232 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2233 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2234 v = a->a; jj = a->j; 2235 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2236 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2237 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2238 } 2239 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2240 PetscFunctionReturn(0); 2241 } 2242 2243 #undef __FUNCT__ 2244 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2245 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2246 { 2247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2248 PetscErrorCode ierr; 2249 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2250 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2251 const PetscInt *irow,*icol; 2252 PetscInt nrows,ncols; 2253 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2254 MatScalar *a_new,*mat_a; 2255 Mat C; 2256 PetscBool stride; 2257 2258 PetscFunctionBegin; 2259 2260 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2261 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2262 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2263 2264 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2265 if (stride) { 2266 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2267 } else { 2268 first = 0; 2269 step = 0; 2270 } 2271 if (stride && step == 1) { 2272 /* special case of contiguous rows */ 2273 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2274 /* loop over new rows determining lens and starting points */ 2275 for (i=0; i<nrows; i++) { 2276 kstart = ai[irow[i]]; 2277 kend = kstart + ailen[irow[i]]; 2278 starts[i] = kstart; 2279 for (k=kstart; k<kend; k++) { 2280 if (aj[k] >= first) { 2281 starts[i] = k; 2282 break; 2283 } 2284 } 2285 sum = 0; 2286 while (k < kend) { 2287 if (aj[k++] >= first+ncols) break; 2288 sum++; 2289 } 2290 lens[i] = sum; 2291 } 2292 /* create submatrix */ 2293 if (scall == MAT_REUSE_MATRIX) { 2294 PetscInt n_cols,n_rows; 2295 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2296 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2297 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2298 C = *B; 2299 } else { 2300 PetscInt rbs,cbs; 2301 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2302 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2303 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2304 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2305 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2306 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2307 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2308 } 2309 c = (Mat_SeqAIJ*)C->data; 2310 2311 /* loop over rows inserting into submatrix */ 2312 a_new = c->a; 2313 j_new = c->j; 2314 i_new = c->i; 2315 2316 for (i=0; i<nrows; i++) { 2317 ii = starts[i]; 2318 lensi = lens[i]; 2319 for (k=0; k<lensi; k++) { 2320 *j_new++ = aj[ii+k] - first; 2321 } 2322 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2323 a_new += lensi; 2324 i_new[i+1] = i_new[i] + lensi; 2325 c->ilen[i] = lensi; 2326 } 2327 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2328 } else { 2329 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2330 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2331 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2332 for (i=0; i<ncols; i++) { 2333 #if defined(PETSC_USE_DEBUG) 2334 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); 2335 #endif 2336 smap[icol[i]] = i+1; 2337 } 2338 2339 /* determine lens of each row */ 2340 for (i=0; i<nrows; i++) { 2341 kstart = ai[irow[i]]; 2342 kend = kstart + a->ilen[irow[i]]; 2343 lens[i] = 0; 2344 for (k=kstart; k<kend; k++) { 2345 if (smap[aj[k]]) { 2346 lens[i]++; 2347 } 2348 } 2349 } 2350 /* Create and fill new matrix */ 2351 if (scall == MAT_REUSE_MATRIX) { 2352 PetscBool equal; 2353 2354 c = (Mat_SeqAIJ*)((*B)->data); 2355 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2356 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2357 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2358 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2359 C = *B; 2360 } else { 2361 PetscInt rbs,cbs; 2362 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2363 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2364 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2365 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2366 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2367 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2368 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2369 } 2370 c = (Mat_SeqAIJ*)(C->data); 2371 for (i=0; i<nrows; i++) { 2372 row = irow[i]; 2373 kstart = ai[row]; 2374 kend = kstart + a->ilen[row]; 2375 mat_i = c->i[i]; 2376 mat_j = c->j + mat_i; 2377 mat_a = c->a + mat_i; 2378 mat_ilen = c->ilen + i; 2379 for (k=kstart; k<kend; k++) { 2380 if ((tcol=smap[a->j[k]])) { 2381 *mat_j++ = tcol - 1; 2382 *mat_a++ = a->a[k]; 2383 (*mat_ilen)++; 2384 2385 } 2386 } 2387 } 2388 /* Free work space */ 2389 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2390 ierr = PetscFree(smap);CHKERRQ(ierr); 2391 ierr = PetscFree(lens);CHKERRQ(ierr); 2392 /* sort */ 2393 for (i = 0; i < nrows; i++) { 2394 PetscInt ilen; 2395 2396 mat_i = c->i[i]; 2397 mat_j = c->j + mat_i; 2398 mat_a = c->a + mat_i; 2399 ilen = c->ilen[i]; 2400 ierr = PetscSortIntWithMatScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2401 } 2402 } 2403 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2404 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2405 2406 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2407 *B = C; 2408 PetscFunctionReturn(0); 2409 } 2410 2411 #undef __FUNCT__ 2412 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2413 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2414 { 2415 PetscErrorCode ierr; 2416 Mat B; 2417 2418 PetscFunctionBegin; 2419 if (scall == MAT_INITIAL_MATRIX) { 2420 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2421 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2422 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2423 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2424 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2425 *subMat = B; 2426 } else { 2427 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2428 } 2429 PetscFunctionReturn(0); 2430 } 2431 2432 #undef __FUNCT__ 2433 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2434 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2435 { 2436 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2437 PetscErrorCode ierr; 2438 Mat outA; 2439 PetscBool row_identity,col_identity; 2440 2441 PetscFunctionBegin; 2442 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2443 2444 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2445 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2446 2447 outA = inA; 2448 outA->factortype = MAT_FACTOR_LU; 2449 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2450 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2451 2452 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2453 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2454 2455 a->row = row; 2456 2457 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2458 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2459 2460 a->col = col; 2461 2462 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2463 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2464 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2465 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2466 2467 if (!a->solve_work) { /* this matrix may have been factored before */ 2468 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2469 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2470 } 2471 2472 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2473 if (row_identity && col_identity) { 2474 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2475 } else { 2476 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2477 } 2478 PetscFunctionReturn(0); 2479 } 2480 2481 #undef __FUNCT__ 2482 #define __FUNCT__ "MatScale_SeqAIJ" 2483 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2484 { 2485 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2486 PetscScalar oalpha = alpha; 2487 PetscErrorCode ierr; 2488 PetscBLASInt one = 1,bnz; 2489 2490 PetscFunctionBegin; 2491 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2492 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2493 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2494 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 #undef __FUNCT__ 2499 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2500 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2501 { 2502 PetscErrorCode ierr; 2503 PetscInt i; 2504 2505 PetscFunctionBegin; 2506 if (scall == MAT_INITIAL_MATRIX) { 2507 ierr = PetscMalloc1(n+1,B);CHKERRQ(ierr); 2508 } 2509 2510 for (i=0; i<n; i++) { 2511 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2512 } 2513 PetscFunctionReturn(0); 2514 } 2515 2516 #undef __FUNCT__ 2517 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2518 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2519 { 2520 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2521 PetscErrorCode ierr; 2522 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2523 const PetscInt *idx; 2524 PetscInt start,end,*ai,*aj; 2525 PetscBT table; 2526 2527 PetscFunctionBegin; 2528 m = A->rmap->n; 2529 ai = a->i; 2530 aj = a->j; 2531 2532 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2533 2534 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2535 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2536 2537 for (i=0; i<is_max; i++) { 2538 /* Initialize the two local arrays */ 2539 isz = 0; 2540 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2541 2542 /* Extract the indices, assume there can be duplicate entries */ 2543 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2544 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2545 2546 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2547 for (j=0; j<n; ++j) { 2548 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2549 } 2550 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2551 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2552 2553 k = 0; 2554 for (j=0; j<ov; j++) { /* for each overlap */ 2555 n = isz; 2556 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2557 row = nidx[k]; 2558 start = ai[row]; 2559 end = ai[row+1]; 2560 for (l = start; l<end; l++) { 2561 val = aj[l]; 2562 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2563 } 2564 } 2565 } 2566 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2567 } 2568 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2569 ierr = PetscFree(nidx);CHKERRQ(ierr); 2570 PetscFunctionReturn(0); 2571 } 2572 2573 /* -------------------------------------------------------------- */ 2574 #undef __FUNCT__ 2575 #define __FUNCT__ "MatPermute_SeqAIJ" 2576 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2577 { 2578 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2579 PetscErrorCode ierr; 2580 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2581 const PetscInt *row,*col; 2582 PetscInt *cnew,j,*lens; 2583 IS icolp,irowp; 2584 PetscInt *cwork = NULL; 2585 PetscScalar *vwork = NULL; 2586 2587 PetscFunctionBegin; 2588 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2589 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2590 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2591 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2592 2593 /* determine lengths of permuted rows */ 2594 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2595 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2596 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2597 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2598 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2599 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2600 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2601 ierr = PetscFree(lens);CHKERRQ(ierr); 2602 2603 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2604 for (i=0; i<m; i++) { 2605 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2606 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2607 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2608 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2609 } 2610 ierr = PetscFree(cnew);CHKERRQ(ierr); 2611 2612 (*B)->assembled = PETSC_FALSE; 2613 2614 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2615 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2616 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2617 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2618 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2619 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 2623 #undef __FUNCT__ 2624 #define __FUNCT__ "MatCopy_SeqAIJ" 2625 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2626 { 2627 PetscErrorCode ierr; 2628 2629 PetscFunctionBegin; 2630 /* If the two matrices have the same copy implementation, use fast copy. */ 2631 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2632 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2633 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2634 2635 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"); 2636 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2637 } else { 2638 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2639 } 2640 PetscFunctionReturn(0); 2641 } 2642 2643 #undef __FUNCT__ 2644 #define __FUNCT__ "MatSetUp_SeqAIJ" 2645 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2646 { 2647 PetscErrorCode ierr; 2648 2649 PetscFunctionBegin; 2650 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2651 PetscFunctionReturn(0); 2652 } 2653 2654 #undef __FUNCT__ 2655 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2656 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2657 { 2658 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2659 2660 PetscFunctionBegin; 2661 *array = a->a; 2662 PetscFunctionReturn(0); 2663 } 2664 2665 #undef __FUNCT__ 2666 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2667 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2668 { 2669 PetscFunctionBegin; 2670 PetscFunctionReturn(0); 2671 } 2672 2673 /* 2674 Computes the number of nonzeros per row needed for preallocation when X and Y 2675 have different nonzero structure. 2676 */ 2677 #undef __FUNCT__ 2678 #define __FUNCT__ "MatAXPYGetPreallocation_SeqX_private" 2679 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2680 { 2681 PetscInt i,j,k,nzx,nzy; 2682 2683 PetscFunctionBegin; 2684 /* Set the number of nonzeros in the new matrix */ 2685 for (i=0; i<m; i++) { 2686 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2687 nzx = xi[i+1] - xi[i]; 2688 nzy = yi[i+1] - yi[i]; 2689 nnz[i] = 0; 2690 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2691 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2692 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2693 nnz[i]++; 2694 } 2695 for (; k<nzy; k++) nnz[i]++; 2696 } 2697 PetscFunctionReturn(0); 2698 } 2699 2700 #undef __FUNCT__ 2701 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2702 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2703 { 2704 PetscInt m = Y->rmap->N; 2705 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2706 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2707 PetscErrorCode ierr; 2708 2709 PetscFunctionBegin; 2710 /* Set the number of nonzeros in the new matrix */ 2711 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2712 PetscFunctionReturn(0); 2713 } 2714 2715 #undef __FUNCT__ 2716 #define __FUNCT__ "MatAXPY_SeqAIJ" 2717 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2718 { 2719 PetscErrorCode ierr; 2720 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2721 PetscBLASInt one=1,bnz; 2722 2723 PetscFunctionBegin; 2724 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2725 if (str == SAME_NONZERO_PATTERN) { 2726 PetscScalar alpha = a; 2727 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2728 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2729 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2730 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2731 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2732 } else { 2733 Mat B; 2734 PetscInt *nnz; 2735 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2736 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2737 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2738 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2739 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2740 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2741 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2742 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2743 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2744 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2745 ierr = PetscFree(nnz);CHKERRQ(ierr); 2746 } 2747 PetscFunctionReturn(0); 2748 } 2749 2750 #undef __FUNCT__ 2751 #define __FUNCT__ "MatConjugate_SeqAIJ" 2752 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2753 { 2754 #if defined(PETSC_USE_COMPLEX) 2755 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2756 PetscInt i,nz; 2757 PetscScalar *a; 2758 2759 PetscFunctionBegin; 2760 nz = aij->nz; 2761 a = aij->a; 2762 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2763 #else 2764 PetscFunctionBegin; 2765 #endif 2766 PetscFunctionReturn(0); 2767 } 2768 2769 #undef __FUNCT__ 2770 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2771 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2772 { 2773 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2774 PetscErrorCode ierr; 2775 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2776 PetscReal atmp; 2777 PetscScalar *x; 2778 MatScalar *aa; 2779 2780 PetscFunctionBegin; 2781 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2782 aa = a->a; 2783 ai = a->i; 2784 aj = a->j; 2785 2786 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2787 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2788 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2789 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2790 for (i=0; i<m; i++) { 2791 ncols = ai[1] - ai[0]; ai++; 2792 x[i] = 0.0; 2793 for (j=0; j<ncols; j++) { 2794 atmp = PetscAbsScalar(*aa); 2795 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2796 aa++; aj++; 2797 } 2798 } 2799 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2800 PetscFunctionReturn(0); 2801 } 2802 2803 #undef __FUNCT__ 2804 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2805 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2806 { 2807 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2808 PetscErrorCode ierr; 2809 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2810 PetscScalar *x; 2811 MatScalar *aa; 2812 2813 PetscFunctionBegin; 2814 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2815 aa = a->a; 2816 ai = a->i; 2817 aj = a->j; 2818 2819 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2820 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2821 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2822 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2823 for (i=0; i<m; i++) { 2824 ncols = ai[1] - ai[0]; ai++; 2825 if (ncols == A->cmap->n) { /* row is dense */ 2826 x[i] = *aa; if (idx) idx[i] = 0; 2827 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2828 x[i] = 0.0; 2829 if (idx) { 2830 idx[i] = 0; /* in case ncols is zero */ 2831 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2832 if (aj[j] > j) { 2833 idx[i] = j; 2834 break; 2835 } 2836 } 2837 } 2838 } 2839 for (j=0; j<ncols; j++) { 2840 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2841 aa++; aj++; 2842 } 2843 } 2844 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2845 PetscFunctionReturn(0); 2846 } 2847 2848 #undef __FUNCT__ 2849 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2850 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2851 { 2852 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2853 PetscErrorCode ierr; 2854 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2855 PetscReal atmp; 2856 PetscScalar *x; 2857 MatScalar *aa; 2858 2859 PetscFunctionBegin; 2860 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2861 aa = a->a; 2862 ai = a->i; 2863 aj = a->j; 2864 2865 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2866 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2867 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2868 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); 2869 for (i=0; i<m; i++) { 2870 ncols = ai[1] - ai[0]; ai++; 2871 if (ncols) { 2872 /* Get first nonzero */ 2873 for (j = 0; j < ncols; j++) { 2874 atmp = PetscAbsScalar(aa[j]); 2875 if (atmp > 1.0e-12) { 2876 x[i] = atmp; 2877 if (idx) idx[i] = aj[j]; 2878 break; 2879 } 2880 } 2881 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2882 } else { 2883 x[i] = 0.0; if (idx) idx[i] = 0; 2884 } 2885 for (j = 0; j < ncols; j++) { 2886 atmp = PetscAbsScalar(*aa); 2887 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2888 aa++; aj++; 2889 } 2890 } 2891 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2892 PetscFunctionReturn(0); 2893 } 2894 2895 #undef __FUNCT__ 2896 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2897 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2898 { 2899 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2900 PetscErrorCode ierr; 2901 PetscInt i,j,m = A->rmap->n,ncols,n; 2902 const PetscInt *ai,*aj; 2903 PetscScalar *x; 2904 const MatScalar *aa; 2905 2906 PetscFunctionBegin; 2907 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2908 aa = a->a; 2909 ai = a->i; 2910 aj = a->j; 2911 2912 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2913 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2914 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2915 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2916 for (i=0; i<m; i++) { 2917 ncols = ai[1] - ai[0]; ai++; 2918 if (ncols == A->cmap->n) { /* row is dense */ 2919 x[i] = *aa; if (idx) idx[i] = 0; 2920 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2921 x[i] = 0.0; 2922 if (idx) { /* find first implicit 0.0 in the row */ 2923 idx[i] = 0; /* in case ncols is zero */ 2924 for (j=0; j<ncols; j++) { 2925 if (aj[j] > j) { 2926 idx[i] = j; 2927 break; 2928 } 2929 } 2930 } 2931 } 2932 for (j=0; j<ncols; j++) { 2933 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2934 aa++; aj++; 2935 } 2936 } 2937 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2938 PetscFunctionReturn(0); 2939 } 2940 2941 #include <petscblaslapack.h> 2942 #include <petsc/private/kernels/blockinvert.h> 2943 2944 #undef __FUNCT__ 2945 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 2946 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2947 { 2948 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2949 PetscErrorCode ierr; 2950 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2951 MatScalar *diag,work[25],*v_work; 2952 PetscReal shift = 0.0; 2953 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2954 2955 PetscFunctionBegin; 2956 allowzeropivot = PetscNot(A->erroriffailure); 2957 if (a->ibdiagvalid) { 2958 if (values) *values = a->ibdiag; 2959 PetscFunctionReturn(0); 2960 } 2961 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2962 if (!a->ibdiag) { 2963 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2964 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2965 } 2966 diag = a->ibdiag; 2967 if (values) *values = a->ibdiag; 2968 /* factor and invert each block */ 2969 switch (bs) { 2970 case 1: 2971 for (i=0; i<mbs; i++) { 2972 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2973 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 2974 if (allowzeropivot) { 2975 A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2976 ierr = PetscInfo1(A,"Zero pivot, row %D\n",i);CHKERRQ(ierr); 2977 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",i); 2978 } 2979 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2980 } 2981 break; 2982 case 2: 2983 for (i=0; i<mbs; i++) { 2984 ij[0] = 2*i; ij[1] = 2*i + 1; 2985 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2986 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2987 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2988 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2989 diag += 4; 2990 } 2991 break; 2992 case 3: 2993 for (i=0; i<mbs; i++) { 2994 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2995 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2996 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2997 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2998 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2999 diag += 9; 3000 } 3001 break; 3002 case 4: 3003 for (i=0; i<mbs; i++) { 3004 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3005 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3006 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3007 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3008 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3009 diag += 16; 3010 } 3011 break; 3012 case 5: 3013 for (i=0; i<mbs; i++) { 3014 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3015 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3016 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3017 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3018 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3019 diag += 25; 3020 } 3021 break; 3022 case 6: 3023 for (i=0; i<mbs; i++) { 3024 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; 3025 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3026 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3027 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3028 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3029 diag += 36; 3030 } 3031 break; 3032 case 7: 3033 for (i=0; i<mbs; i++) { 3034 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; 3035 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3036 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3037 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3038 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3039 diag += 49; 3040 } 3041 break; 3042 default: 3043 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3044 for (i=0; i<mbs; i++) { 3045 for (j=0; j<bs; j++) { 3046 IJ[j] = bs*i + j; 3047 } 3048 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3049 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3050 if (zeropivotdetected) A->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3051 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3052 diag += bs2; 3053 } 3054 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3055 } 3056 a->ibdiagvalid = PETSC_TRUE; 3057 PetscFunctionReturn(0); 3058 } 3059 3060 #undef __FUNCT__ 3061 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3062 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3063 { 3064 PetscErrorCode ierr; 3065 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3066 PetscScalar a; 3067 PetscInt m,n,i,j,col; 3068 3069 PetscFunctionBegin; 3070 if (!x->assembled) { 3071 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3072 for (i=0; i<m; i++) { 3073 for (j=0; j<aij->imax[i]; j++) { 3074 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3075 col = (PetscInt)(n*PetscRealPart(a)); 3076 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3077 } 3078 } 3079 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3080 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3081 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3082 PetscFunctionReturn(0); 3083 } 3084 3085 #undef __FUNCT__ 3086 #define __FUNCT__ "MatShift_SeqAIJ" 3087 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3088 { 3089 PetscErrorCode ierr; 3090 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3091 3092 PetscFunctionBegin; 3093 if (!Y->preallocated || !aij->nz) { 3094 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3095 } 3096 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3097 PetscFunctionReturn(0); 3098 } 3099 3100 /* -------------------------------------------------------------------*/ 3101 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3102 MatGetRow_SeqAIJ, 3103 MatRestoreRow_SeqAIJ, 3104 MatMult_SeqAIJ, 3105 /* 4*/ MatMultAdd_SeqAIJ, 3106 MatMultTranspose_SeqAIJ, 3107 MatMultTransposeAdd_SeqAIJ, 3108 0, 3109 0, 3110 0, 3111 /* 10*/ 0, 3112 MatLUFactor_SeqAIJ, 3113 0, 3114 MatSOR_SeqAIJ, 3115 MatTranspose_SeqAIJ, 3116 /*1 5*/ MatGetInfo_SeqAIJ, 3117 MatEqual_SeqAIJ, 3118 MatGetDiagonal_SeqAIJ, 3119 MatDiagonalScale_SeqAIJ, 3120 MatNorm_SeqAIJ, 3121 /* 20*/ 0, 3122 MatAssemblyEnd_SeqAIJ, 3123 MatSetOption_SeqAIJ, 3124 MatZeroEntries_SeqAIJ, 3125 /* 24*/ MatZeroRows_SeqAIJ, 3126 0, 3127 0, 3128 0, 3129 0, 3130 /* 29*/ MatSetUp_SeqAIJ, 3131 0, 3132 0, 3133 0, 3134 0, 3135 /* 34*/ MatDuplicate_SeqAIJ, 3136 0, 3137 0, 3138 MatILUFactor_SeqAIJ, 3139 0, 3140 /* 39*/ MatAXPY_SeqAIJ, 3141 MatGetSubMatrices_SeqAIJ, 3142 MatIncreaseOverlap_SeqAIJ, 3143 MatGetValues_SeqAIJ, 3144 MatCopy_SeqAIJ, 3145 /* 44*/ MatGetRowMax_SeqAIJ, 3146 MatScale_SeqAIJ, 3147 MatShift_SeqAIJ, 3148 MatDiagonalSet_SeqAIJ, 3149 MatZeroRowsColumns_SeqAIJ, 3150 /* 49*/ MatSetRandom_SeqAIJ, 3151 MatGetRowIJ_SeqAIJ, 3152 MatRestoreRowIJ_SeqAIJ, 3153 MatGetColumnIJ_SeqAIJ, 3154 MatRestoreColumnIJ_SeqAIJ, 3155 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3156 0, 3157 0, 3158 MatPermute_SeqAIJ, 3159 0, 3160 /* 59*/ 0, 3161 MatDestroy_SeqAIJ, 3162 MatView_SeqAIJ, 3163 0, 3164 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3165 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3166 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3167 0, 3168 0, 3169 0, 3170 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3171 MatGetRowMinAbs_SeqAIJ, 3172 0, 3173 MatSetColoring_SeqAIJ, 3174 0, 3175 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3176 MatFDColoringApply_AIJ, 3177 0, 3178 0, 3179 0, 3180 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3181 0, 3182 0, 3183 0, 3184 MatLoad_SeqAIJ, 3185 /* 84*/ MatIsSymmetric_SeqAIJ, 3186 MatIsHermitian_SeqAIJ, 3187 0, 3188 0, 3189 0, 3190 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3191 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3192 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3193 MatPtAP_SeqAIJ_SeqAIJ, 3194 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3195 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3196 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3197 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3198 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3199 0, 3200 /* 99*/ 0, 3201 0, 3202 0, 3203 MatConjugate_SeqAIJ, 3204 0, 3205 /*104*/ MatSetValuesRow_SeqAIJ, 3206 MatRealPart_SeqAIJ, 3207 MatImaginaryPart_SeqAIJ, 3208 0, 3209 0, 3210 /*109*/ MatMatSolve_SeqAIJ, 3211 0, 3212 MatGetRowMin_SeqAIJ, 3213 0, 3214 MatMissingDiagonal_SeqAIJ, 3215 /*114*/ 0, 3216 0, 3217 0, 3218 0, 3219 0, 3220 /*119*/ 0, 3221 0, 3222 0, 3223 0, 3224 MatGetMultiProcBlock_SeqAIJ, 3225 /*124*/ MatFindNonzeroRows_SeqAIJ, 3226 MatGetColumnNorms_SeqAIJ, 3227 MatInvertBlockDiagonal_SeqAIJ, 3228 0, 3229 0, 3230 /*129*/ 0, 3231 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3232 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3233 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3234 MatTransposeColoringCreate_SeqAIJ, 3235 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3236 MatTransColoringApplyDenToSp_SeqAIJ, 3237 MatRARt_SeqAIJ_SeqAIJ, 3238 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3239 MatRARtNumeric_SeqAIJ_SeqAIJ, 3240 /*139*/0, 3241 0, 3242 0, 3243 MatFDColoringSetUp_SeqXAIJ, 3244 MatFindOffBlockDiagonalEntries_SeqAIJ, 3245 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ 3246 }; 3247 3248 #undef __FUNCT__ 3249 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3250 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3251 { 3252 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3253 PetscInt i,nz,n; 3254 3255 PetscFunctionBegin; 3256 nz = aij->maxnz; 3257 n = mat->rmap->n; 3258 for (i=0; i<nz; i++) { 3259 aij->j[i] = indices[i]; 3260 } 3261 aij->nz = nz; 3262 for (i=0; i<n; i++) { 3263 aij->ilen[i] = aij->imax[i]; 3264 } 3265 PetscFunctionReturn(0); 3266 } 3267 3268 #undef __FUNCT__ 3269 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3270 /*@ 3271 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3272 in the matrix. 3273 3274 Input Parameters: 3275 + mat - the SeqAIJ matrix 3276 - indices - the column indices 3277 3278 Level: advanced 3279 3280 Notes: 3281 This can be called if you have precomputed the nonzero structure of the 3282 matrix and want to provide it to the matrix object to improve the performance 3283 of the MatSetValues() operation. 3284 3285 You MUST have set the correct numbers of nonzeros per row in the call to 3286 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3287 3288 MUST be called before any calls to MatSetValues(); 3289 3290 The indices should start with zero, not one. 3291 3292 @*/ 3293 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3294 { 3295 PetscErrorCode ierr; 3296 3297 PetscFunctionBegin; 3298 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3299 PetscValidPointer(indices,2); 3300 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3301 PetscFunctionReturn(0); 3302 } 3303 3304 /* ----------------------------------------------------------------------------------------*/ 3305 3306 #undef __FUNCT__ 3307 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3308 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3309 { 3310 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3311 PetscErrorCode ierr; 3312 size_t nz = aij->i[mat->rmap->n]; 3313 3314 PetscFunctionBegin; 3315 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3316 3317 /* allocate space for values if not already there */ 3318 if (!aij->saved_values) { 3319 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3320 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3321 } 3322 3323 /* copy values over */ 3324 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3325 PetscFunctionReturn(0); 3326 } 3327 3328 #undef __FUNCT__ 3329 #define __FUNCT__ "MatStoreValues" 3330 /*@ 3331 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3332 example, reuse of the linear part of a Jacobian, while recomputing the 3333 nonlinear portion. 3334 3335 Collect on Mat 3336 3337 Input Parameters: 3338 . mat - the matrix (currently only AIJ matrices support this option) 3339 3340 Level: advanced 3341 3342 Common Usage, with SNESSolve(): 3343 $ Create Jacobian matrix 3344 $ Set linear terms into matrix 3345 $ Apply boundary conditions to matrix, at this time matrix must have 3346 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3347 $ boundary conditions again will not change the nonzero structure 3348 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3349 $ ierr = MatStoreValues(mat); 3350 $ Call SNESSetJacobian() with matrix 3351 $ In your Jacobian routine 3352 $ ierr = MatRetrieveValues(mat); 3353 $ Set nonlinear terms in matrix 3354 3355 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3356 $ // build linear portion of Jacobian 3357 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3358 $ ierr = MatStoreValues(mat); 3359 $ loop over nonlinear iterations 3360 $ ierr = MatRetrieveValues(mat); 3361 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3362 $ // call MatAssemblyBegin/End() on matrix 3363 $ Solve linear system with Jacobian 3364 $ endloop 3365 3366 Notes: 3367 Matrix must already be assemblied before calling this routine 3368 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3369 calling this routine. 3370 3371 When this is called multiple times it overwrites the previous set of stored values 3372 and does not allocated additional space. 3373 3374 .seealso: MatRetrieveValues() 3375 3376 @*/ 3377 PetscErrorCode MatStoreValues(Mat mat) 3378 { 3379 PetscErrorCode ierr; 3380 3381 PetscFunctionBegin; 3382 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3383 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3384 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3385 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3386 PetscFunctionReturn(0); 3387 } 3388 3389 #undef __FUNCT__ 3390 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3391 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3392 { 3393 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3394 PetscErrorCode ierr; 3395 PetscInt nz = aij->i[mat->rmap->n]; 3396 3397 PetscFunctionBegin; 3398 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3399 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3400 /* copy values over */ 3401 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3402 PetscFunctionReturn(0); 3403 } 3404 3405 #undef __FUNCT__ 3406 #define __FUNCT__ "MatRetrieveValues" 3407 /*@ 3408 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3409 example, reuse of the linear part of a Jacobian, while recomputing the 3410 nonlinear portion. 3411 3412 Collect on Mat 3413 3414 Input Parameters: 3415 . mat - the matrix (currently on AIJ matrices support this option) 3416 3417 Level: advanced 3418 3419 .seealso: MatStoreValues() 3420 3421 @*/ 3422 PetscErrorCode MatRetrieveValues(Mat mat) 3423 { 3424 PetscErrorCode ierr; 3425 3426 PetscFunctionBegin; 3427 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3428 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3429 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3430 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3431 PetscFunctionReturn(0); 3432 } 3433 3434 3435 /* --------------------------------------------------------------------------------*/ 3436 #undef __FUNCT__ 3437 #define __FUNCT__ "MatCreateSeqAIJ" 3438 /*@C 3439 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3440 (the default parallel PETSc format). For good matrix assembly performance 3441 the user should preallocate the matrix storage by setting the parameter nz 3442 (or the array nnz). By setting these parameters accurately, performance 3443 during matrix assembly can be increased by more than a factor of 50. 3444 3445 Collective on MPI_Comm 3446 3447 Input Parameters: 3448 + comm - MPI communicator, set to PETSC_COMM_SELF 3449 . m - number of rows 3450 . n - number of columns 3451 . nz - number of nonzeros per row (same for all rows) 3452 - nnz - array containing the number of nonzeros in the various rows 3453 (possibly different for each row) or NULL 3454 3455 Output Parameter: 3456 . A - the matrix 3457 3458 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3459 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3460 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3461 3462 Notes: 3463 If nnz is given then nz is ignored 3464 3465 The AIJ format (also called the Yale sparse matrix format or 3466 compressed row storage), is fully compatible with standard Fortran 77 3467 storage. That is, the stored row and column indices can begin at 3468 either one (as in Fortran) or zero. See the users' manual for details. 3469 3470 Specify the preallocated storage with either nz or nnz (not both). 3471 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3472 allocation. For large problems you MUST preallocate memory or you 3473 will get TERRIBLE performance, see the users' manual chapter on matrices. 3474 3475 By default, this format uses inodes (identical nodes) when possible, to 3476 improve numerical efficiency of matrix-vector products and solves. We 3477 search for consecutive rows with the same nonzero structure, thereby 3478 reusing matrix information to achieve increased efficiency. 3479 3480 Options Database Keys: 3481 + -mat_no_inode - Do not use inodes 3482 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3483 3484 Level: intermediate 3485 3486 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3487 3488 @*/ 3489 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3490 { 3491 PetscErrorCode ierr; 3492 3493 PetscFunctionBegin; 3494 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3495 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3496 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3497 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3498 PetscFunctionReturn(0); 3499 } 3500 3501 #undef __FUNCT__ 3502 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3503 /*@C 3504 MatSeqAIJSetPreallocation - For good matrix assembly performance 3505 the user should preallocate the matrix storage by setting the parameter nz 3506 (or the array nnz). By setting these parameters accurately, performance 3507 during matrix assembly can be increased by more than a factor of 50. 3508 3509 Collective on MPI_Comm 3510 3511 Input Parameters: 3512 + B - The matrix 3513 . nz - number of nonzeros per row (same for all rows) 3514 - nnz - array containing the number of nonzeros in the various rows 3515 (possibly different for each row) or NULL 3516 3517 Notes: 3518 If nnz is given then nz is ignored 3519 3520 The AIJ format (also called the Yale sparse matrix format or 3521 compressed row storage), is fully compatible with standard Fortran 77 3522 storage. That is, the stored row and column indices can begin at 3523 either one (as in Fortran) or zero. See the users' manual for details. 3524 3525 Specify the preallocated storage with either nz or nnz (not both). 3526 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3527 allocation. For large problems you MUST preallocate memory or you 3528 will get TERRIBLE performance, see the users' manual chapter on matrices. 3529 3530 You can call MatGetInfo() to get information on how effective the preallocation was; 3531 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3532 You can also run with the option -info and look for messages with the string 3533 malloc in them to see if additional memory allocation was needed. 3534 3535 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3536 entries or columns indices 3537 3538 By default, this format uses inodes (identical nodes) when possible, to 3539 improve numerical efficiency of matrix-vector products and solves. We 3540 search for consecutive rows with the same nonzero structure, thereby 3541 reusing matrix information to achieve increased efficiency. 3542 3543 Options Database Keys: 3544 + -mat_no_inode - Do not use inodes 3545 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3546 - -mat_aij_oneindex - Internally use indexing starting at 1 3547 rather than 0. Note that when calling MatSetValues(), 3548 the user still MUST index entries starting at 0! 3549 3550 Level: intermediate 3551 3552 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3553 3554 @*/ 3555 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3556 { 3557 PetscErrorCode ierr; 3558 3559 PetscFunctionBegin; 3560 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3561 PetscValidType(B,1); 3562 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3563 PetscFunctionReturn(0); 3564 } 3565 3566 #undef __FUNCT__ 3567 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3568 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3569 { 3570 Mat_SeqAIJ *b; 3571 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3572 PetscErrorCode ierr; 3573 PetscInt i; 3574 3575 PetscFunctionBegin; 3576 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3577 if (nz == MAT_SKIP_ALLOCATION) { 3578 skipallocation = PETSC_TRUE; 3579 nz = 0; 3580 } 3581 3582 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3583 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3584 3585 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3586 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3587 if (nnz) { 3588 for (i=0; i<B->rmap->n; i++) { 3589 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]); 3590 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); 3591 } 3592 } 3593 3594 B->preallocated = PETSC_TRUE; 3595 3596 b = (Mat_SeqAIJ*)B->data; 3597 3598 if (!skipallocation) { 3599 if (!b->imax) { 3600 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3601 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3602 } 3603 if (!nnz) { 3604 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3605 else if (nz < 0) nz = 1; 3606 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3607 nz = nz*B->rmap->n; 3608 } else { 3609 nz = 0; 3610 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3611 } 3612 /* b->ilen will count nonzeros in each row so far. */ 3613 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3614 3615 /* allocate the matrix space */ 3616 /* FIXME: should B's old memory be unlogged? */ 3617 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3618 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3619 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3620 b->i[0] = 0; 3621 for (i=1; i<B->rmap->n+1; i++) { 3622 b->i[i] = b->i[i-1] + b->imax[i-1]; 3623 } 3624 b->singlemalloc = PETSC_TRUE; 3625 b->free_a = PETSC_TRUE; 3626 b->free_ij = PETSC_TRUE; 3627 } else { 3628 b->free_a = PETSC_FALSE; 3629 b->free_ij = PETSC_FALSE; 3630 } 3631 3632 b->nz = 0; 3633 b->maxnz = nz; 3634 B->info.nz_unneeded = (double)b->maxnz; 3635 if (realalloc) { 3636 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3637 } 3638 PetscFunctionReturn(0); 3639 } 3640 3641 #undef __FUNCT__ 3642 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3643 /*@ 3644 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3645 3646 Input Parameters: 3647 + B - the matrix 3648 . i - the indices into j for the start of each row (starts with zero) 3649 . j - the column indices for each row (starts with zero) these must be sorted for each row 3650 - v - optional values in the matrix 3651 3652 Level: developer 3653 3654 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3655 3656 .keywords: matrix, aij, compressed row, sparse, sequential 3657 3658 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3659 @*/ 3660 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3661 { 3662 PetscErrorCode ierr; 3663 3664 PetscFunctionBegin; 3665 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3666 PetscValidType(B,1); 3667 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3668 PetscFunctionReturn(0); 3669 } 3670 3671 #undef __FUNCT__ 3672 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3673 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3674 { 3675 PetscInt i; 3676 PetscInt m,n; 3677 PetscInt nz; 3678 PetscInt *nnz, nz_max = 0; 3679 PetscScalar *values; 3680 PetscErrorCode ierr; 3681 3682 PetscFunctionBegin; 3683 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3684 3685 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3686 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3687 3688 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3689 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3690 for (i = 0; i < m; i++) { 3691 nz = Ii[i+1]- Ii[i]; 3692 nz_max = PetscMax(nz_max, nz); 3693 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3694 nnz[i] = nz; 3695 } 3696 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3697 ierr = PetscFree(nnz);CHKERRQ(ierr); 3698 3699 if (v) { 3700 values = (PetscScalar*) v; 3701 } else { 3702 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3703 } 3704 3705 for (i = 0; i < m; i++) { 3706 nz = Ii[i+1] - Ii[i]; 3707 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3708 } 3709 3710 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3711 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3712 3713 if (!v) { 3714 ierr = PetscFree(values);CHKERRQ(ierr); 3715 } 3716 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3717 PetscFunctionReturn(0); 3718 } 3719 3720 #include <../src/mat/impls/dense/seq/dense.h> 3721 #include <petsc/private/kernels/petscaxpy.h> 3722 3723 #undef __FUNCT__ 3724 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3725 /* 3726 Computes (B'*A')' since computing B*A directly is untenable 3727 3728 n p p 3729 ( ) ( ) ( ) 3730 m ( A ) * n ( B ) = m ( C ) 3731 ( ) ( ) ( ) 3732 3733 */ 3734 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3735 { 3736 PetscErrorCode ierr; 3737 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3738 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3739 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3740 PetscInt i,n,m,q,p; 3741 const PetscInt *ii,*idx; 3742 const PetscScalar *b,*a,*a_q; 3743 PetscScalar *c,*c_q; 3744 3745 PetscFunctionBegin; 3746 m = A->rmap->n; 3747 n = A->cmap->n; 3748 p = B->cmap->n; 3749 a = sub_a->v; 3750 b = sub_b->a; 3751 c = sub_c->v; 3752 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3753 3754 ii = sub_b->i; 3755 idx = sub_b->j; 3756 for (i=0; i<n; i++) { 3757 q = ii[i+1] - ii[i]; 3758 while (q-->0) { 3759 c_q = c + m*(*idx); 3760 a_q = a + m*i; 3761 PetscKernelAXPY(c_q,*b,a_q,m); 3762 idx++; 3763 b++; 3764 } 3765 } 3766 PetscFunctionReturn(0); 3767 } 3768 3769 #undef __FUNCT__ 3770 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3771 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3772 { 3773 PetscErrorCode ierr; 3774 PetscInt m=A->rmap->n,n=B->cmap->n; 3775 Mat Cmat; 3776 3777 PetscFunctionBegin; 3778 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); 3779 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3780 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3781 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3782 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3783 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3784 3785 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3786 3787 *C = Cmat; 3788 PetscFunctionReturn(0); 3789 } 3790 3791 /* ----------------------------------------------------------------*/ 3792 #undef __FUNCT__ 3793 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3794 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3795 { 3796 PetscErrorCode ierr; 3797 3798 PetscFunctionBegin; 3799 if (scall == MAT_INITIAL_MATRIX) { 3800 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3801 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3802 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3803 } 3804 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3805 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3806 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3807 PetscFunctionReturn(0); 3808 } 3809 3810 3811 /*MC 3812 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3813 based on compressed sparse row format. 3814 3815 Options Database Keys: 3816 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3817 3818 Level: beginner 3819 3820 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3821 M*/ 3822 3823 /*MC 3824 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3825 3826 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3827 and MATMPIAIJ otherwise. As a result, for single process communicators, 3828 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3829 for communicators controlling multiple processes. It is recommended that you call both of 3830 the above preallocation routines for simplicity. 3831 3832 Options Database Keys: 3833 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3834 3835 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3836 enough exist. 3837 3838 Level: beginner 3839 3840 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3841 M*/ 3842 3843 /*MC 3844 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3845 3846 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3847 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3848 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3849 for communicators controlling multiple processes. It is recommended that you call both of 3850 the above preallocation routines for simplicity. 3851 3852 Options Database Keys: 3853 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3854 3855 Level: beginner 3856 3857 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3858 M*/ 3859 3860 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3861 #if defined(PETSC_HAVE_ELEMENTAL) 3862 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3863 #endif 3864 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3865 3866 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3867 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3868 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3869 #endif 3870 3871 3872 #undef __FUNCT__ 3873 #define __FUNCT__ "MatSeqAIJGetArray" 3874 /*@C 3875 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3876 3877 Not Collective 3878 3879 Input Parameter: 3880 . mat - a MATSEQAIJ matrix 3881 3882 Output Parameter: 3883 . array - pointer to the data 3884 3885 Level: intermediate 3886 3887 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3888 @*/ 3889 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3890 { 3891 PetscErrorCode ierr; 3892 3893 PetscFunctionBegin; 3894 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3895 PetscFunctionReturn(0); 3896 } 3897 3898 #undef __FUNCT__ 3899 #define __FUNCT__ "MatSeqAIJGetMaxRowNonzeros" 3900 /*@C 3901 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3902 3903 Not Collective 3904 3905 Input Parameter: 3906 . mat - a MATSEQAIJ matrix 3907 3908 Output Parameter: 3909 . nz - the maximum number of nonzeros in any row 3910 3911 Level: intermediate 3912 3913 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3914 @*/ 3915 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3916 { 3917 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3918 3919 PetscFunctionBegin; 3920 *nz = aij->rmax; 3921 PetscFunctionReturn(0); 3922 } 3923 3924 #undef __FUNCT__ 3925 #define __FUNCT__ "MatSeqAIJRestoreArray" 3926 /*@C 3927 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3928 3929 Not Collective 3930 3931 Input Parameters: 3932 . mat - a MATSEQAIJ matrix 3933 . array - pointer to the data 3934 3935 Level: intermediate 3936 3937 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3938 @*/ 3939 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3940 { 3941 PetscErrorCode ierr; 3942 3943 PetscFunctionBegin; 3944 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3945 PetscFunctionReturn(0); 3946 } 3947 3948 #undef __FUNCT__ 3949 #define __FUNCT__ "MatCreate_SeqAIJ" 3950 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3951 { 3952 Mat_SeqAIJ *b; 3953 PetscErrorCode ierr; 3954 PetscMPIInt size; 3955 3956 PetscFunctionBegin; 3957 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3958 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3959 3960 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3961 3962 B->data = (void*)b; 3963 3964 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3965 3966 b->row = 0; 3967 b->col = 0; 3968 b->icol = 0; 3969 b->reallocs = 0; 3970 b->ignorezeroentries = PETSC_FALSE; 3971 b->roworiented = PETSC_TRUE; 3972 b->nonew = 0; 3973 b->diag = 0; 3974 b->solve_work = 0; 3975 B->spptr = 0; 3976 b->saved_values = 0; 3977 b->idiag = 0; 3978 b->mdiag = 0; 3979 b->ssor_work = 0; 3980 b->omega = 1.0; 3981 b->fshift = 0.0; 3982 b->idiagvalid = PETSC_FALSE; 3983 b->ibdiagvalid = PETSC_FALSE; 3984 b->keepnonzeropattern = PETSC_FALSE; 3985 3986 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3987 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3988 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3989 3990 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3991 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3992 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3993 #endif 3994 3995 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3996 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3997 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3998 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3999 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4000 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4001 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4002 #if defined(PETSC_HAVE_ELEMENTAL) 4003 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4004 #endif 4005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4009 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4011 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4012 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4013 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4014 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4015 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4016 PetscFunctionReturn(0); 4017 } 4018 4019 #undef __FUNCT__ 4020 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4021 /* 4022 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4023 */ 4024 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4025 { 4026 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4027 PetscErrorCode ierr; 4028 PetscInt i,m = A->rmap->n; 4029 4030 PetscFunctionBegin; 4031 c = (Mat_SeqAIJ*)C->data; 4032 4033 C->factortype = A->factortype; 4034 c->row = 0; 4035 c->col = 0; 4036 c->icol = 0; 4037 c->reallocs = 0; 4038 4039 C->assembled = PETSC_TRUE; 4040 4041 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4042 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4043 4044 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4045 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4046 for (i=0; i<m; i++) { 4047 c->imax[i] = a->imax[i]; 4048 c->ilen[i] = a->ilen[i]; 4049 } 4050 4051 /* allocate the matrix space */ 4052 if (mallocmatspace) { 4053 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4054 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4055 4056 c->singlemalloc = PETSC_TRUE; 4057 4058 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4059 if (m > 0) { 4060 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4061 if (cpvalues == MAT_COPY_VALUES) { 4062 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4063 } else { 4064 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4065 } 4066 } 4067 } 4068 4069 c->ignorezeroentries = a->ignorezeroentries; 4070 c->roworiented = a->roworiented; 4071 c->nonew = a->nonew; 4072 if (a->diag) { 4073 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4074 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4075 for (i=0; i<m; i++) { 4076 c->diag[i] = a->diag[i]; 4077 } 4078 } else c->diag = 0; 4079 4080 c->solve_work = 0; 4081 c->saved_values = 0; 4082 c->idiag = 0; 4083 c->ssor_work = 0; 4084 c->keepnonzeropattern = a->keepnonzeropattern; 4085 c->free_a = PETSC_TRUE; 4086 c->free_ij = PETSC_TRUE; 4087 4088 c->rmax = a->rmax; 4089 c->nz = a->nz; 4090 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4091 C->preallocated = PETSC_TRUE; 4092 4093 c->compressedrow.use = a->compressedrow.use; 4094 c->compressedrow.nrows = a->compressedrow.nrows; 4095 if (a->compressedrow.use) { 4096 i = a->compressedrow.nrows; 4097 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4098 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4099 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4100 } else { 4101 c->compressedrow.use = PETSC_FALSE; 4102 c->compressedrow.i = NULL; 4103 c->compressedrow.rindex = NULL; 4104 } 4105 c->nonzerorowcnt = a->nonzerorowcnt; 4106 C->nonzerostate = A->nonzerostate; 4107 4108 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4109 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4110 PetscFunctionReturn(0); 4111 } 4112 4113 #undef __FUNCT__ 4114 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4115 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4116 { 4117 PetscErrorCode ierr; 4118 4119 PetscFunctionBegin; 4120 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4121 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4122 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4123 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4124 } 4125 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4126 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4127 PetscFunctionReturn(0); 4128 } 4129 4130 #undef __FUNCT__ 4131 #define __FUNCT__ "MatLoad_SeqAIJ" 4132 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4133 { 4134 Mat_SeqAIJ *a; 4135 PetscErrorCode ierr; 4136 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4137 int fd; 4138 PetscMPIInt size; 4139 MPI_Comm comm; 4140 PetscInt bs = newMat->rmap->bs; 4141 4142 PetscFunctionBegin; 4143 /* force binary viewer to load .info file if it has not yet done so */ 4144 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4145 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4146 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4147 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4148 4149 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4150 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4151 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4152 if (bs < 0) bs = 1; 4153 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4154 4155 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4156 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4157 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4158 M = header[1]; N = header[2]; nz = header[3]; 4159 4160 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4161 4162 /* read in row lengths */ 4163 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4164 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4165 4166 /* check if sum of rowlengths is same as nz */ 4167 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4168 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); 4169 4170 /* set global size if not set already*/ 4171 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4172 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4173 } else { 4174 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4175 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4176 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4177 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4178 } 4179 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); 4180 } 4181 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4182 a = (Mat_SeqAIJ*)newMat->data; 4183 4184 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4185 4186 /* read in nonzero values */ 4187 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4188 4189 /* set matrix "i" values */ 4190 a->i[0] = 0; 4191 for (i=1; i<= M; i++) { 4192 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4193 a->ilen[i-1] = rowlengths[i-1]; 4194 } 4195 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4196 4197 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4198 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4199 PetscFunctionReturn(0); 4200 } 4201 4202 #undef __FUNCT__ 4203 #define __FUNCT__ "MatEqual_SeqAIJ" 4204 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4205 { 4206 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4207 PetscErrorCode ierr; 4208 #if defined(PETSC_USE_COMPLEX) 4209 PetscInt k; 4210 #endif 4211 4212 PetscFunctionBegin; 4213 /* If the matrix dimensions are not equal,or no of nonzeros */ 4214 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4215 *flg = PETSC_FALSE; 4216 PetscFunctionReturn(0); 4217 } 4218 4219 /* if the a->i are the same */ 4220 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4221 if (!*flg) PetscFunctionReturn(0); 4222 4223 /* if a->j are the same */ 4224 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4225 if (!*flg) PetscFunctionReturn(0); 4226 4227 /* if a->a are the same */ 4228 #if defined(PETSC_USE_COMPLEX) 4229 for (k=0; k<a->nz; k++) { 4230 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4231 *flg = PETSC_FALSE; 4232 PetscFunctionReturn(0); 4233 } 4234 } 4235 #else 4236 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4237 #endif 4238 PetscFunctionReturn(0); 4239 } 4240 4241 #undef __FUNCT__ 4242 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4243 /*@ 4244 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4245 provided by the user. 4246 4247 Collective on MPI_Comm 4248 4249 Input Parameters: 4250 + comm - must be an MPI communicator of size 1 4251 . m - number of rows 4252 . n - number of columns 4253 . i - row indices 4254 . j - column indices 4255 - a - matrix values 4256 4257 Output Parameter: 4258 . mat - the matrix 4259 4260 Level: intermediate 4261 4262 Notes: 4263 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4264 once the matrix is destroyed and not before 4265 4266 You cannot set new nonzero locations into this matrix, that will generate an error. 4267 4268 The i and j indices are 0 based 4269 4270 The format which is used for the sparse matrix input, is equivalent to a 4271 row-major ordering.. i.e for the following matrix, the input data expected is 4272 as shown 4273 4274 $ 1 0 0 4275 $ 2 0 3 4276 $ 4 5 6 4277 $ 4278 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4279 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4280 $ v = {1,2,3,4,5,6} [size = 6] 4281 4282 4283 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4284 4285 @*/ 4286 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4287 { 4288 PetscErrorCode ierr; 4289 PetscInt ii; 4290 Mat_SeqAIJ *aij; 4291 #if defined(PETSC_USE_DEBUG) 4292 PetscInt jj; 4293 #endif 4294 4295 PetscFunctionBegin; 4296 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4297 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4298 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4299 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4300 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4301 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4302 aij = (Mat_SeqAIJ*)(*mat)->data; 4303 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4304 4305 aij->i = i; 4306 aij->j = j; 4307 aij->a = a; 4308 aij->singlemalloc = PETSC_FALSE; 4309 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4310 aij->free_a = PETSC_FALSE; 4311 aij->free_ij = PETSC_FALSE; 4312 4313 for (ii=0; ii<m; ii++) { 4314 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4315 #if defined(PETSC_USE_DEBUG) 4316 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]); 4317 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4318 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); 4319 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); 4320 } 4321 #endif 4322 } 4323 #if defined(PETSC_USE_DEBUG) 4324 for (ii=0; ii<aij->i[m]; ii++) { 4325 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4326 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]); 4327 } 4328 #endif 4329 4330 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4331 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4332 PetscFunctionReturn(0); 4333 } 4334 #undef __FUNCT__ 4335 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4336 /*@C 4337 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4338 provided by the user. 4339 4340 Collective on MPI_Comm 4341 4342 Input Parameters: 4343 + comm - must be an MPI communicator of size 1 4344 . m - number of rows 4345 . n - number of columns 4346 . i - row indices 4347 . j - column indices 4348 . a - matrix values 4349 . nz - number of nonzeros 4350 - idx - 0 or 1 based 4351 4352 Output Parameter: 4353 . mat - the matrix 4354 4355 Level: intermediate 4356 4357 Notes: 4358 The i and j indices are 0 based 4359 4360 The format which is used for the sparse matrix input, is equivalent to a 4361 row-major ordering.. i.e for the following matrix, the input data expected is 4362 as shown: 4363 4364 1 0 0 4365 2 0 3 4366 4 5 6 4367 4368 i = {0,1,1,2,2,2} 4369 j = {0,0,2,0,1,2} 4370 v = {1,2,3,4,5,6} 4371 4372 4373 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4374 4375 @*/ 4376 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4377 { 4378 PetscErrorCode ierr; 4379 PetscInt ii, *nnz, one = 1,row,col; 4380 4381 4382 PetscFunctionBegin; 4383 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4384 for (ii = 0; ii < nz; ii++) { 4385 nnz[i[ii] - !!idx] += 1; 4386 } 4387 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4388 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4389 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4390 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4391 for (ii = 0; ii < nz; ii++) { 4392 if (idx) { 4393 row = i[ii] - 1; 4394 col = j[ii] - 1; 4395 } else { 4396 row = i[ii]; 4397 col = j[ii]; 4398 } 4399 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4400 } 4401 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4402 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4403 ierr = PetscFree(nnz);CHKERRQ(ierr); 4404 PetscFunctionReturn(0); 4405 } 4406 4407 #undef __FUNCT__ 4408 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4409 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4410 { 4411 PetscErrorCode ierr; 4412 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4413 4414 PetscFunctionBegin; 4415 if (coloring->ctype == IS_COLORING_GLOBAL) { 4416 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4417 a->coloring = coloring; 4418 } else if (coloring->ctype == IS_COLORING_LOCAL) { 4419 PetscInt i,*larray; 4420 ISColoring ocoloring; 4421 ISColoringValue *colors; 4422 4423 /* set coloring for diagonal portion */ 4424 ierr = PetscMalloc1(A->cmap->n,&larray);CHKERRQ(ierr); 4425 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4426 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4427 ierr = PetscMalloc1(A->cmap->n,&colors);CHKERRQ(ierr); 4428 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4429 ierr = PetscFree(larray);CHKERRQ(ierr); 4430 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,PETSC_OWN_POINTER,&ocoloring);CHKERRQ(ierr); 4431 a->coloring = ocoloring; 4432 } 4433 PetscFunctionReturn(0); 4434 } 4435 4436 #undef __FUNCT__ 4437 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4438 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4439 { 4440 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4441 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4442 MatScalar *v = a->a; 4443 PetscScalar *values = (PetscScalar*)advalues; 4444 ISColoringValue *color; 4445 4446 PetscFunctionBegin; 4447 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4448 color = a->coloring->colors; 4449 /* loop over rows */ 4450 for (i=0; i<m; i++) { 4451 nz = ii[i+1] - ii[i]; 4452 /* loop over columns putting computed value into matrix */ 4453 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4454 values += nl; /* jump to next row of derivatives */ 4455 } 4456 PetscFunctionReturn(0); 4457 } 4458 4459 #undef __FUNCT__ 4460 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4461 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4462 { 4463 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4464 PetscErrorCode ierr; 4465 4466 PetscFunctionBegin; 4467 a->idiagvalid = PETSC_FALSE; 4468 a->ibdiagvalid = PETSC_FALSE; 4469 4470 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4471 PetscFunctionReturn(0); 4472 } 4473 4474 #undef __FUNCT__ 4475 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqAIJ" 4476 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4477 { 4478 PetscErrorCode ierr; 4479 4480 PetscFunctionBegin; 4481 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4482 PetscFunctionReturn(0); 4483 } 4484 4485 /* 4486 Permute A into C's *local* index space using rowemb,colemb. 4487 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4488 of [0,m), colemb is in [0,n). 4489 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4490 */ 4491 #undef __FUNCT__ 4492 #define __FUNCT__ "MatSetSeqMat_SeqAIJ" 4493 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4494 { 4495 /* If making this function public, change the error returned in this function away from _PLIB. */ 4496 PetscErrorCode ierr; 4497 Mat_SeqAIJ *Baij; 4498 PetscBool seqaij; 4499 PetscInt m,n,*nz,i,j,count; 4500 PetscScalar v; 4501 const PetscInt *rowindices,*colindices; 4502 4503 PetscFunctionBegin; 4504 if (!B) PetscFunctionReturn(0); 4505 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4506 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4507 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4508 if (rowemb) { 4509 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4510 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); 4511 } else { 4512 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4513 } 4514 if (colemb) { 4515 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4516 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); 4517 } else { 4518 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4519 } 4520 4521 Baij = (Mat_SeqAIJ*)(B->data); 4522 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4523 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4524 for (i=0; i<B->rmap->n; i++) { 4525 nz[i] = Baij->i[i+1] - Baij->i[i]; 4526 } 4527 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4528 ierr = PetscFree(nz);CHKERRQ(ierr); 4529 } 4530 if (pattern == SUBSET_NONZERO_PATTERN) { 4531 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4532 } 4533 count = 0; 4534 rowindices = NULL; 4535 colindices = NULL; 4536 if (rowemb) { 4537 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4538 } 4539 if (colemb) { 4540 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4541 } 4542 for (i=0; i<B->rmap->n; i++) { 4543 PetscInt row; 4544 row = i; 4545 if (rowindices) row = rowindices[i]; 4546 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4547 PetscInt col; 4548 col = Baij->j[count]; 4549 if (colindices) col = colindices[col]; 4550 v = Baij->a[count]; 4551 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4552 ++count; 4553 } 4554 } 4555 /* FIXME: set C's nonzerostate correctly. */ 4556 /* Assembly for C is necessary. */ 4557 C->preallocated = PETSC_TRUE; 4558 C->assembled = PETSC_TRUE; 4559 C->was_assembled = PETSC_FALSE; 4560 PetscFunctionReturn(0); 4561 } 4562 4563 4564 /* 4565 Special version for direct calls from Fortran 4566 */ 4567 #include <petsc/private/fortranimpl.h> 4568 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4569 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4570 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4571 #define matsetvaluesseqaij_ matsetvaluesseqaij 4572 #endif 4573 4574 /* Change these macros so can be used in void function */ 4575 #undef CHKERRQ 4576 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4577 #undef SETERRQ2 4578 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4579 #undef SETERRQ3 4580 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4581 4582 #undef __FUNCT__ 4583 #define __FUNCT__ "matsetvaluesseqaij_" 4584 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) 4585 { 4586 Mat A = *AA; 4587 PetscInt m = *mm, n = *nn; 4588 InsertMode is = *isis; 4589 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4590 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4591 PetscInt *imax,*ai,*ailen; 4592 PetscErrorCode ierr; 4593 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4594 MatScalar *ap,value,*aa; 4595 PetscBool ignorezeroentries = a->ignorezeroentries; 4596 PetscBool roworiented = a->roworiented; 4597 4598 PetscFunctionBegin; 4599 MatCheckPreallocated(A,1); 4600 imax = a->imax; 4601 ai = a->i; 4602 ailen = a->ilen; 4603 aj = a->j; 4604 aa = a->a; 4605 4606 for (k=0; k<m; k++) { /* loop over added rows */ 4607 row = im[k]; 4608 if (row < 0) continue; 4609 #if defined(PETSC_USE_DEBUG) 4610 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4611 #endif 4612 rp = aj + ai[row]; ap = aa + ai[row]; 4613 rmax = imax[row]; nrow = ailen[row]; 4614 low = 0; 4615 high = nrow; 4616 for (l=0; l<n; l++) { /* loop over added columns */ 4617 if (in[l] < 0) continue; 4618 #if defined(PETSC_USE_DEBUG) 4619 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4620 #endif 4621 col = in[l]; 4622 if (roworiented) value = v[l + k*n]; 4623 else value = v[k + l*m]; 4624 4625 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4626 4627 if (col <= lastcol) low = 0; 4628 else high = nrow; 4629 lastcol = col; 4630 while (high-low > 5) { 4631 t = (low+high)/2; 4632 if (rp[t] > col) high = t; 4633 else low = t; 4634 } 4635 for (i=low; i<high; i++) { 4636 if (rp[i] > col) break; 4637 if (rp[i] == col) { 4638 if (is == ADD_VALUES) ap[i] += value; 4639 else ap[i] = value; 4640 goto noinsert; 4641 } 4642 } 4643 if (value == 0.0 && ignorezeroentries) goto noinsert; 4644 if (nonew == 1) goto noinsert; 4645 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4646 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4647 N = nrow++ - 1; a->nz++; high++; 4648 /* shift up all the later entries in this row */ 4649 for (ii=N; ii>=i; ii--) { 4650 rp[ii+1] = rp[ii]; 4651 ap[ii+1] = ap[ii]; 4652 } 4653 rp[i] = col; 4654 ap[i] = value; 4655 A->nonzerostate++; 4656 noinsert:; 4657 low = i + 1; 4658 } 4659 ailen[row] = nrow; 4660 } 4661 PetscFunctionReturnVoid(); 4662 } 4663 4664