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