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