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