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