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