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