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 *y,*z; 1157 const PetscScalar *x; 1158 const MatScalar *aa; 1159 PetscErrorCode ierr; 1160 PetscInt m = A->rmap->n,*aj,*ii; 1161 PetscInt n,i,*ridx=PETSC_NULL; 1162 PetscScalar sum; 1163 PetscBool usecprow=a->compressedrow.use; 1164 1165 PetscFunctionBegin; 1166 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1167 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1168 if (zz != yy) { 1169 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1170 } else { 1171 z = y; 1172 } 1173 1174 aj = a->j; 1175 aa = a->a; 1176 ii = a->i; 1177 if (usecprow){ /* use compressed row format */ 1178 if (zz != yy){ 1179 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1180 } 1181 m = a->compressedrow.nrows; 1182 ii = a->compressedrow.i; 1183 ridx = a->compressedrow.rindex; 1184 for (i=0; i<m; i++){ 1185 n = ii[i+1] - ii[i]; 1186 aj = a->j + ii[i]; 1187 aa = a->a + ii[i]; 1188 sum = y[*ridx]; 1189 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1190 z[*ridx++] = sum; 1191 } 1192 } else { /* do not use compressed row format */ 1193 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1194 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1195 #else 1196 for (i=0; i<m; i++) { 1197 n = ii[i+1] - ii[i]; 1198 aj = a->j + ii[i]; 1199 aa = a->a + ii[i]; 1200 sum = y[i]; 1201 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1202 z[i] = sum; 1203 } 1204 #endif 1205 } 1206 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1207 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1208 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1209 if (zz != yy) { 1210 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1211 } 1212 #if defined(PETSC_HAVE_CUSP) 1213 /* 1214 ierr = VecView(xx,0);CHKERRQ(ierr); 1215 ierr = VecView(zz,0);CHKERRQ(ierr); 1216 ierr = MatView(A,0);CHKERRQ(ierr); 1217 */ 1218 #endif 1219 PetscFunctionReturn(0); 1220 } 1221 1222 /* 1223 Adds diagonal pointers to sparse matrix structure. 1224 */ 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1227 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1228 { 1229 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1230 PetscErrorCode ierr; 1231 PetscInt i,j,m = A->rmap->n; 1232 1233 PetscFunctionBegin; 1234 if (!a->diag) { 1235 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1236 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1237 } 1238 for (i=0; i<A->rmap->n; i++) { 1239 a->diag[i] = a->i[i+1]; 1240 for (j=a->i[i]; j<a->i[i+1]; j++) { 1241 if (a->j[j] == i) { 1242 a->diag[i] = j; 1243 break; 1244 } 1245 } 1246 } 1247 PetscFunctionReturn(0); 1248 } 1249 1250 /* 1251 Checks for missing diagonals 1252 */ 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1255 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1256 { 1257 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1258 PetscInt *diag,*jj = a->j,i; 1259 1260 PetscFunctionBegin; 1261 *missing = PETSC_FALSE; 1262 if (A->rmap->n > 0 && !jj) { 1263 *missing = PETSC_TRUE; 1264 if (d) *d = 0; 1265 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1266 } else { 1267 diag = a->diag; 1268 for (i=0; i<A->rmap->n; i++) { 1269 if (jj[diag[i]] != i) { 1270 *missing = PETSC_TRUE; 1271 if (d) *d = i; 1272 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1273 } 1274 } 1275 } 1276 PetscFunctionReturn(0); 1277 } 1278 1279 EXTERN_C_BEGIN 1280 #undef __FUNCT__ 1281 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1282 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1283 { 1284 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1285 PetscErrorCode ierr; 1286 PetscInt i,*diag,m = A->rmap->n; 1287 MatScalar *v = a->a; 1288 PetscScalar *idiag,*mdiag; 1289 1290 PetscFunctionBegin; 1291 if (a->idiagvalid) PetscFunctionReturn(0); 1292 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1293 diag = a->diag; 1294 if (!a->idiag) { 1295 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1296 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1297 v = a->a; 1298 } 1299 mdiag = a->mdiag; 1300 idiag = a->idiag; 1301 1302 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1303 for (i=0; i<m; i++) { 1304 mdiag[i] = v[diag[i]]; 1305 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1306 idiag[i] = 1.0/v[diag[i]]; 1307 } 1308 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1309 } else { 1310 for (i=0; i<m; i++) { 1311 mdiag[i] = v[diag[i]]; 1312 idiag[i] = omega/(fshift + v[diag[i]]); 1313 } 1314 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1315 } 1316 a->idiagvalid = PETSC_TRUE; 1317 PetscFunctionReturn(0); 1318 } 1319 EXTERN_C_END 1320 1321 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1322 #undef __FUNCT__ 1323 #define __FUNCT__ "MatSOR_SeqAIJ" 1324 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1325 { 1326 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1327 PetscScalar *x,d,sum,*t,scale; 1328 const MatScalar *v = a->a,*idiag=0,*mdiag; 1329 const PetscScalar *b, *bs,*xb, *ts; 1330 PetscErrorCode ierr; 1331 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1332 const PetscInt *idx,*diag; 1333 1334 PetscFunctionBegin; 1335 its = its*lits; 1336 1337 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1338 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1339 a->fshift = fshift; 1340 a->omega = omega; 1341 1342 diag = a->diag; 1343 t = a->ssor_work; 1344 idiag = a->idiag; 1345 mdiag = a->mdiag; 1346 1347 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1348 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1349 CHKMEMQ; 1350 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1351 if (flag == SOR_APPLY_UPPER) { 1352 /* apply (U + D/omega) to the vector */ 1353 bs = b; 1354 for (i=0; i<m; i++) { 1355 d = fshift + mdiag[i]; 1356 n = a->i[i+1] - diag[i] - 1; 1357 idx = a->j + diag[i] + 1; 1358 v = a->a + diag[i] + 1; 1359 sum = b[i]*d/omega; 1360 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1361 x[i] = sum; 1362 } 1363 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1364 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1365 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1366 PetscFunctionReturn(0); 1367 } 1368 1369 if (flag == SOR_APPLY_LOWER) { 1370 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1371 } else if (flag & SOR_EISENSTAT) { 1372 /* Let A = L + U + D; where L is lower trianglar, 1373 U is upper triangular, E = D/omega; This routine applies 1374 1375 (L + E)^{-1} A (U + E)^{-1} 1376 1377 to a vector efficiently using Eisenstat's trick. 1378 */ 1379 scale = (2.0/omega) - 1.0; 1380 1381 /* x = (E + U)^{-1} b */ 1382 for (i=m-1; i>=0; i--) { 1383 n = a->i[i+1] - diag[i] - 1; 1384 idx = a->j + diag[i] + 1; 1385 v = a->a + diag[i] + 1; 1386 sum = b[i]; 1387 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1388 x[i] = sum*idiag[i]; 1389 } 1390 1391 /* t = b - (2*E - D)x */ 1392 v = a->a; 1393 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1394 1395 /* t = (E + L)^{-1}t */ 1396 ts = t; 1397 diag = a->diag; 1398 for (i=0; i<m; i++) { 1399 n = diag[i] - a->i[i]; 1400 idx = a->j + a->i[i]; 1401 v = a->a + a->i[i]; 1402 sum = t[i]; 1403 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1404 t[i] = sum*idiag[i]; 1405 /* x = x + t */ 1406 x[i] += t[i]; 1407 } 1408 1409 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1410 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1411 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 if (flag & SOR_ZERO_INITIAL_GUESS) { 1415 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1416 for (i=0; i<m; i++) { 1417 n = diag[i] - a->i[i]; 1418 idx = a->j + a->i[i]; 1419 v = a->a + a->i[i]; 1420 sum = b[i]; 1421 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1422 t[i] = sum; 1423 x[i] = sum*idiag[i]; 1424 } 1425 xb = t; 1426 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1427 } else xb = b; 1428 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1429 for (i=m-1; i>=0; i--) { 1430 n = a->i[i+1] - diag[i] - 1; 1431 idx = a->j + diag[i] + 1; 1432 v = a->a + diag[i] + 1; 1433 sum = xb[i]; 1434 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1435 if (xb == b) { 1436 x[i] = sum*idiag[i]; 1437 } else { 1438 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1439 } 1440 } 1441 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1442 } 1443 its--; 1444 } 1445 while (its--) { 1446 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1447 for (i=0; i<m; i++) { 1448 n = a->i[i+1] - a->i[i]; 1449 idx = a->j + a->i[i]; 1450 v = a->a + a->i[i]; 1451 sum = b[i]; 1452 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1453 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1454 } 1455 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1456 } 1457 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1458 for (i=m-1; i>=0; i--) { 1459 n = a->i[i+1] - a->i[i]; 1460 idx = a->j + a->i[i]; 1461 v = a->a + a->i[i]; 1462 sum = b[i]; 1463 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1464 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1465 } 1466 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1467 } 1468 } 1469 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1470 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1471 CHKMEMQ; PetscFunctionReturn(0); 1472 } 1473 1474 1475 #undef __FUNCT__ 1476 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1477 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1478 { 1479 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1480 1481 PetscFunctionBegin; 1482 info->block_size = 1.0; 1483 info->nz_allocated = (double)a->maxnz; 1484 info->nz_used = (double)a->nz; 1485 info->nz_unneeded = (double)(a->maxnz - a->nz); 1486 info->assemblies = (double)A->num_ass; 1487 info->mallocs = (double)A->info.mallocs; 1488 info->memory = ((PetscObject)A)->mem; 1489 if (A->factortype) { 1490 info->fill_ratio_given = A->info.fill_ratio_given; 1491 info->fill_ratio_needed = A->info.fill_ratio_needed; 1492 info->factor_mallocs = A->info.factor_mallocs; 1493 } else { 1494 info->fill_ratio_given = 0; 1495 info->fill_ratio_needed = 0; 1496 info->factor_mallocs = 0; 1497 } 1498 PetscFunctionReturn(0); 1499 } 1500 1501 #undef __FUNCT__ 1502 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1503 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1504 { 1505 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1506 PetscInt i,m = A->rmap->n - 1,d = 0; 1507 PetscErrorCode ierr; 1508 const PetscScalar *xx; 1509 PetscScalar *bb; 1510 PetscBool missing; 1511 1512 PetscFunctionBegin; 1513 if (x && b) { 1514 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1515 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1516 for (i=0; i<N; i++) { 1517 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1518 bb[rows[i]] = diag*xx[rows[i]]; 1519 } 1520 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1521 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1522 } 1523 1524 if (a->keepnonzeropattern) { 1525 for (i=0; i<N; i++) { 1526 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1527 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1528 } 1529 if (diag != 0.0) { 1530 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1531 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1532 for (i=0; i<N; i++) { 1533 a->a[a->diag[rows[i]]] = diag; 1534 } 1535 } 1536 A->same_nonzero = PETSC_TRUE; 1537 } else { 1538 if (diag != 0.0) { 1539 for (i=0; i<N; i++) { 1540 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1541 if (a->ilen[rows[i]] > 0) { 1542 a->ilen[rows[i]] = 1; 1543 a->a[a->i[rows[i]]] = diag; 1544 a->j[a->i[rows[i]]] = rows[i]; 1545 } else { /* in case row was completely empty */ 1546 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1547 } 1548 } 1549 } else { 1550 for (i=0; i<N; i++) { 1551 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1552 a->ilen[rows[i]] = 0; 1553 } 1554 } 1555 A->same_nonzero = PETSC_FALSE; 1556 } 1557 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1563 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1564 { 1565 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1566 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1567 PetscErrorCode ierr; 1568 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1569 const PetscScalar *xx; 1570 PetscScalar *bb; 1571 1572 PetscFunctionBegin; 1573 if (x && b) { 1574 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1575 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1576 vecs = PETSC_TRUE; 1577 } 1578 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1579 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1580 for (i=0; i<N; i++) { 1581 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1582 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1583 zeroed[rows[i]] = PETSC_TRUE; 1584 } 1585 for (i=0; i<A->rmap->n; i++) { 1586 if (!zeroed[i]) { 1587 for (j=a->i[i]; j<a->i[i+1]; j++) { 1588 if (zeroed[a->j[j]]) { 1589 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1590 a->a[j] = 0.0; 1591 } 1592 } 1593 } else if (vecs) bb[i] = diag*xx[i]; 1594 } 1595 if (x && b) { 1596 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1597 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1598 } 1599 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1600 if (diag != 0.0) { 1601 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1602 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1603 for (i=0; i<N; i++) { 1604 a->a[a->diag[rows[i]]] = diag; 1605 } 1606 } 1607 A->same_nonzero = PETSC_TRUE; 1608 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1609 PetscFunctionReturn(0); 1610 } 1611 1612 #undef __FUNCT__ 1613 #define __FUNCT__ "MatGetRow_SeqAIJ" 1614 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1615 { 1616 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1617 PetscInt *itmp; 1618 1619 PetscFunctionBegin; 1620 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1621 1622 *nz = a->i[row+1] - a->i[row]; 1623 if (v) *v = a->a + a->i[row]; 1624 if (idx) { 1625 itmp = a->j + a->i[row]; 1626 if (*nz) { 1627 *idx = itmp; 1628 } 1629 else *idx = 0; 1630 } 1631 PetscFunctionReturn(0); 1632 } 1633 1634 /* remove this function? */ 1635 #undef __FUNCT__ 1636 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1637 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1638 { 1639 PetscFunctionBegin; 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "MatNorm_SeqAIJ" 1645 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1646 { 1647 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1648 MatScalar *v = a->a; 1649 PetscReal sum = 0.0; 1650 PetscErrorCode ierr; 1651 PetscInt i,j; 1652 1653 PetscFunctionBegin; 1654 if (type == NORM_FROBENIUS) { 1655 for (i=0; i<a->nz; i++) { 1656 #if defined(PETSC_USE_COMPLEX) 1657 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1658 #else 1659 sum += (*v)*(*v); v++; 1660 #endif 1661 } 1662 *nrm = sqrt(sum); 1663 } else if (type == NORM_1) { 1664 PetscReal *tmp; 1665 PetscInt *jj = a->j; 1666 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1667 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1668 *nrm = 0.0; 1669 for (j=0; j<a->nz; j++) { 1670 tmp[*jj++] += PetscAbsScalar(*v); v++; 1671 } 1672 for (j=0; j<A->cmap->n; j++) { 1673 if (tmp[j] > *nrm) *nrm = tmp[j]; 1674 } 1675 ierr = PetscFree(tmp);CHKERRQ(ierr); 1676 } else if (type == NORM_INFINITY) { 1677 *nrm = 0.0; 1678 for (j=0; j<A->rmap->n; j++) { 1679 v = a->a + a->i[j]; 1680 sum = 0.0; 1681 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1682 sum += PetscAbsScalar(*v); v++; 1683 } 1684 if (sum > *nrm) *nrm = sum; 1685 } 1686 } else { 1687 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1688 } 1689 PetscFunctionReturn(0); 1690 } 1691 1692 #undef __FUNCT__ 1693 #define __FUNCT__ "MatTranspose_SeqAIJ" 1694 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1695 { 1696 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1697 Mat C; 1698 PetscErrorCode ierr; 1699 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1700 MatScalar *array = a->a; 1701 1702 PetscFunctionBegin; 1703 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"); 1704 1705 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1706 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1707 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1708 1709 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1710 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1711 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1712 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1713 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1714 ierr = PetscFree(col);CHKERRQ(ierr); 1715 } else { 1716 C = *B; 1717 } 1718 1719 for (i=0; i<m; i++) { 1720 len = ai[i+1]-ai[i]; 1721 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1722 array += len; 1723 aj += len; 1724 } 1725 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1726 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1727 1728 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1729 *B = C; 1730 } else { 1731 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 1732 } 1733 PetscFunctionReturn(0); 1734 } 1735 1736 EXTERN_C_BEGIN 1737 #undef __FUNCT__ 1738 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1739 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 1740 { 1741 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1742 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1743 MatScalar *va,*vb; 1744 PetscErrorCode ierr; 1745 PetscInt ma,na,mb,nb, i; 1746 1747 PetscFunctionBegin; 1748 bij = (Mat_SeqAIJ *) B->data; 1749 1750 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1751 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1752 if (ma!=nb || na!=mb){ 1753 *f = PETSC_FALSE; 1754 PetscFunctionReturn(0); 1755 } 1756 aii = aij->i; bii = bij->i; 1757 adx = aij->j; bdx = bij->j; 1758 va = aij->a; vb = bij->a; 1759 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1760 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1761 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1762 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1763 1764 *f = PETSC_TRUE; 1765 for (i=0; i<ma; i++) { 1766 while (aptr[i]<aii[i+1]) { 1767 PetscInt idc,idr; 1768 PetscScalar vc,vr; 1769 /* column/row index/value */ 1770 idc = adx[aptr[i]]; 1771 idr = bdx[bptr[idc]]; 1772 vc = va[aptr[i]]; 1773 vr = vb[bptr[idc]]; 1774 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1775 *f = PETSC_FALSE; 1776 goto done; 1777 } else { 1778 aptr[i]++; 1779 if (B || i!=idc) bptr[idc]++; 1780 } 1781 } 1782 } 1783 done: 1784 ierr = PetscFree(aptr);CHKERRQ(ierr); 1785 ierr = PetscFree(bptr);CHKERRQ(ierr); 1786 PetscFunctionReturn(0); 1787 } 1788 EXTERN_C_END 1789 1790 EXTERN_C_BEGIN 1791 #undef __FUNCT__ 1792 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1793 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 1794 { 1795 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1796 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1797 MatScalar *va,*vb; 1798 PetscErrorCode ierr; 1799 PetscInt ma,na,mb,nb, i; 1800 1801 PetscFunctionBegin; 1802 bij = (Mat_SeqAIJ *) B->data; 1803 1804 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1805 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1806 if (ma!=nb || na!=mb){ 1807 *f = PETSC_FALSE; 1808 PetscFunctionReturn(0); 1809 } 1810 aii = aij->i; bii = bij->i; 1811 adx = aij->j; bdx = bij->j; 1812 va = aij->a; vb = bij->a; 1813 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1814 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1815 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1816 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1817 1818 *f = PETSC_TRUE; 1819 for (i=0; i<ma; i++) { 1820 while (aptr[i]<aii[i+1]) { 1821 PetscInt idc,idr; 1822 PetscScalar vc,vr; 1823 /* column/row index/value */ 1824 idc = adx[aptr[i]]; 1825 idr = bdx[bptr[idc]]; 1826 vc = va[aptr[i]]; 1827 vr = vb[bptr[idc]]; 1828 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1829 *f = PETSC_FALSE; 1830 goto done; 1831 } else { 1832 aptr[i]++; 1833 if (B || i!=idc) bptr[idc]++; 1834 } 1835 } 1836 } 1837 done: 1838 ierr = PetscFree(aptr);CHKERRQ(ierr); 1839 ierr = PetscFree(bptr);CHKERRQ(ierr); 1840 PetscFunctionReturn(0); 1841 } 1842 EXTERN_C_END 1843 1844 #undef __FUNCT__ 1845 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1846 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 1847 { 1848 PetscErrorCode ierr; 1849 PetscFunctionBegin; 1850 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1851 PetscFunctionReturn(0); 1852 } 1853 1854 #undef __FUNCT__ 1855 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1856 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 1857 { 1858 PetscErrorCode ierr; 1859 PetscFunctionBegin; 1860 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1861 PetscFunctionReturn(0); 1862 } 1863 1864 #undef __FUNCT__ 1865 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1866 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1867 { 1868 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1869 PetscScalar *l,*r,x; 1870 MatScalar *v; 1871 PetscErrorCode ierr; 1872 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1873 1874 PetscFunctionBegin; 1875 if (ll) { 1876 /* The local size is used so that VecMPI can be passed to this routine 1877 by MatDiagonalScale_MPIAIJ */ 1878 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1879 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1880 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1881 v = a->a; 1882 for (i=0; i<m; i++) { 1883 x = l[i]; 1884 M = a->i[i+1] - a->i[i]; 1885 for (j=0; j<M; j++) { (*v++) *= x;} 1886 } 1887 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1888 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1889 } 1890 if (rr) { 1891 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1892 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1893 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1894 v = a->a; jj = a->j; 1895 for (i=0; i<nz; i++) { 1896 (*v++) *= r[*jj++]; 1897 } 1898 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1899 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1900 } 1901 PetscFunctionReturn(0); 1902 } 1903 1904 #undef __FUNCT__ 1905 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1906 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1907 { 1908 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1909 PetscErrorCode ierr; 1910 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1911 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1912 const PetscInt *irow,*icol; 1913 PetscInt nrows,ncols; 1914 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1915 MatScalar *a_new,*mat_a; 1916 Mat C; 1917 PetscBool stride,sorted; 1918 1919 PetscFunctionBegin; 1920 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1921 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1922 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1923 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1924 1925 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1926 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1927 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1928 1929 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1930 ierr = PetscTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 1931 if (stride && step == 1) { 1932 /* special case of contiguous rows */ 1933 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1934 /* loop over new rows determining lens and starting points */ 1935 for (i=0; i<nrows; i++) { 1936 kstart = ai[irow[i]]; 1937 kend = kstart + ailen[irow[i]]; 1938 for (k=kstart; k<kend; k++) { 1939 if (aj[k] >= first) { 1940 starts[i] = k; 1941 break; 1942 } 1943 } 1944 sum = 0; 1945 while (k < kend) { 1946 if (aj[k++] >= first+ncols) break; 1947 sum++; 1948 } 1949 lens[i] = sum; 1950 } 1951 /* create submatrix */ 1952 if (scall == MAT_REUSE_MATRIX) { 1953 PetscInt n_cols,n_rows; 1954 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1955 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1956 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1957 C = *B; 1958 } else { 1959 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1960 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1961 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1962 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1963 } 1964 c = (Mat_SeqAIJ*)C->data; 1965 1966 /* loop over rows inserting into submatrix */ 1967 a_new = c->a; 1968 j_new = c->j; 1969 i_new = c->i; 1970 1971 for (i=0; i<nrows; i++) { 1972 ii = starts[i]; 1973 lensi = lens[i]; 1974 for (k=0; k<lensi; k++) { 1975 *j_new++ = aj[ii+k] - first; 1976 } 1977 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1978 a_new += lensi; 1979 i_new[i+1] = i_new[i] + lensi; 1980 c->ilen[i] = lensi; 1981 } 1982 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1983 } else { 1984 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1985 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1986 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1987 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1988 for (i=0; i<ncols; i++) { 1989 #if defined(PETSC_USE_DEBUG) 1990 if (icol[i] >= oldcols) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%D] %D <= A->cmap->n %D",i,icol[i],oldcols); 1991 #endif 1992 smap[icol[i]] = i+1; 1993 } 1994 1995 /* determine lens of each row */ 1996 for (i=0; i<nrows; i++) { 1997 kstart = ai[irow[i]]; 1998 kend = kstart + a->ilen[irow[i]]; 1999 lens[i] = 0; 2000 for (k=kstart; k<kend; k++) { 2001 if (smap[aj[k]]) { 2002 lens[i]++; 2003 } 2004 } 2005 } 2006 /* Create and fill new matrix */ 2007 if (scall == MAT_REUSE_MATRIX) { 2008 PetscBool equal; 2009 2010 c = (Mat_SeqAIJ *)((*B)->data); 2011 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2012 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2013 if (!equal) { 2014 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2015 } 2016 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2017 C = *B; 2018 } else { 2019 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2020 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2021 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2022 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2023 } 2024 c = (Mat_SeqAIJ *)(C->data); 2025 for (i=0; i<nrows; i++) { 2026 row = irow[i]; 2027 kstart = ai[row]; 2028 kend = kstart + a->ilen[row]; 2029 mat_i = c->i[i]; 2030 mat_j = c->j + mat_i; 2031 mat_a = c->a + mat_i; 2032 mat_ilen = c->ilen + i; 2033 for (k=kstart; k<kend; k++) { 2034 if ((tcol=smap[a->j[k]])) { 2035 *mat_j++ = tcol - 1; 2036 *mat_a++ = a->a[k]; 2037 (*mat_ilen)++; 2038 2039 } 2040 } 2041 } 2042 /* Free work space */ 2043 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2044 ierr = PetscFree(smap);CHKERRQ(ierr); 2045 ierr = PetscFree(lens);CHKERRQ(ierr); 2046 } 2047 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2048 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2049 2050 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2051 *B = C; 2052 PetscFunctionReturn(0); 2053 } 2054 2055 #undef __FUNCT__ 2056 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2057 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,Mat* subMat) 2058 { 2059 PetscErrorCode ierr; 2060 Mat B; 2061 2062 PetscFunctionBegin; 2063 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2064 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2065 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2066 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2067 *subMat = B; 2068 PetscFunctionReturn(0); 2069 } 2070 2071 #undef __FUNCT__ 2072 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2073 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2074 { 2075 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2076 PetscErrorCode ierr; 2077 Mat outA; 2078 PetscBool row_identity,col_identity; 2079 2080 PetscFunctionBegin; 2081 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2082 2083 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2084 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2085 2086 outA = inA; 2087 outA->factortype = MAT_FACTOR_LU; 2088 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2089 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2090 a->row = row; 2091 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2092 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2093 a->col = col; 2094 2095 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2096 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2097 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2098 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 2099 2100 if (!a->solve_work) { /* this matrix may have been factored before */ 2101 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 2102 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2103 } 2104 2105 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2106 if (row_identity && col_identity) { 2107 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2108 } else { 2109 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2110 } 2111 PetscFunctionReturn(0); 2112 } 2113 2114 #undef __FUNCT__ 2115 #define __FUNCT__ "MatScale_SeqAIJ" 2116 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2117 { 2118 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2119 PetscScalar oalpha = alpha; 2120 PetscErrorCode ierr; 2121 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 2122 2123 PetscFunctionBegin; 2124 BLASscal_(&bnz,&oalpha,a->a,&one); 2125 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2126 PetscFunctionReturn(0); 2127 } 2128 2129 #undef __FUNCT__ 2130 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2131 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2132 { 2133 PetscErrorCode ierr; 2134 PetscInt i; 2135 2136 PetscFunctionBegin; 2137 if (scall == MAT_INITIAL_MATRIX) { 2138 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 2139 } 2140 2141 for (i=0; i<n; i++) { 2142 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2143 } 2144 PetscFunctionReturn(0); 2145 } 2146 2147 #undef __FUNCT__ 2148 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2149 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2150 { 2151 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2152 PetscErrorCode ierr; 2153 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2154 const PetscInt *idx; 2155 PetscInt start,end,*ai,*aj; 2156 PetscBT table; 2157 2158 PetscFunctionBegin; 2159 m = A->rmap->n; 2160 ai = a->i; 2161 aj = a->j; 2162 2163 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2164 2165 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 2166 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 2167 2168 for (i=0; i<is_max; i++) { 2169 /* Initialize the two local arrays */ 2170 isz = 0; 2171 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2172 2173 /* Extract the indices, assume there can be duplicate entries */ 2174 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2175 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2176 2177 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2178 for (j=0; j<n ; ++j){ 2179 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 2180 } 2181 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2182 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2183 2184 k = 0; 2185 for (j=0; j<ov; j++){ /* for each overlap */ 2186 n = isz; 2187 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2188 row = nidx[k]; 2189 start = ai[row]; 2190 end = ai[row+1]; 2191 for (l = start; l<end ; l++){ 2192 val = aj[l] ; 2193 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2194 } 2195 } 2196 } 2197 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2198 } 2199 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2200 ierr = PetscFree(nidx);CHKERRQ(ierr); 2201 PetscFunctionReturn(0); 2202 } 2203 2204 /* -------------------------------------------------------------- */ 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "MatPermute_SeqAIJ" 2207 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2208 { 2209 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2210 PetscErrorCode ierr; 2211 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2212 const PetscInt *row,*col; 2213 PetscInt *cnew,j,*lens; 2214 IS icolp,irowp; 2215 PetscInt *cwork = PETSC_NULL; 2216 PetscScalar *vwork = PETSC_NULL; 2217 2218 PetscFunctionBegin; 2219 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2220 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2221 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2222 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2223 2224 /* determine lengths of permuted rows */ 2225 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2226 for (i=0; i<m; i++) { 2227 lens[row[i]] = a->i[i+1] - a->i[i]; 2228 } 2229 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2230 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2231 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2232 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2233 ierr = PetscFree(lens);CHKERRQ(ierr); 2234 2235 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2236 for (i=0; i<m; i++) { 2237 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2238 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2239 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2240 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2241 } 2242 ierr = PetscFree(cnew);CHKERRQ(ierr); 2243 (*B)->assembled = PETSC_FALSE; 2244 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2245 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2246 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2247 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2248 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2249 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2250 PetscFunctionReturn(0); 2251 } 2252 2253 #undef __FUNCT__ 2254 #define __FUNCT__ "MatCopy_SeqAIJ" 2255 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2256 { 2257 PetscErrorCode ierr; 2258 2259 PetscFunctionBegin; 2260 /* If the two matrices have the same copy implementation, use fast copy. */ 2261 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2262 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2263 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2264 2265 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"); 2266 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2267 } else { 2268 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2269 } 2270 PetscFunctionReturn(0); 2271 } 2272 2273 #undef __FUNCT__ 2274 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2275 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2276 { 2277 PetscErrorCode ierr; 2278 2279 PetscFunctionBegin; 2280 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2281 PetscFunctionReturn(0); 2282 } 2283 2284 #undef __FUNCT__ 2285 #define __FUNCT__ "MatGetArray_SeqAIJ" 2286 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2287 { 2288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2289 PetscFunctionBegin; 2290 *array = a->a; 2291 PetscFunctionReturn(0); 2292 } 2293 2294 #undef __FUNCT__ 2295 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2296 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2297 { 2298 PetscFunctionBegin; 2299 PetscFunctionReturn(0); 2300 } 2301 2302 #undef __FUNCT__ 2303 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2304 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2305 { 2306 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2307 PetscErrorCode ierr; 2308 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2309 PetscScalar dx,*y,*xx,*w3_array; 2310 PetscScalar *vscale_array; 2311 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2312 Vec w1,w2,w3; 2313 void *fctx = coloring->fctx; 2314 PetscBool flg = PETSC_FALSE; 2315 2316 PetscFunctionBegin; 2317 if (!coloring->w1) { 2318 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2319 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2320 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2321 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2322 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2323 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2324 } 2325 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2326 2327 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2328 ierr = PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2329 if (flg) { 2330 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2331 } else { 2332 PetscBool assembled; 2333 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2334 if (assembled) { 2335 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2336 } 2337 } 2338 2339 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2340 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2341 2342 /* 2343 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2344 coloring->F for the coarser grids from the finest 2345 */ 2346 if (coloring->F) { 2347 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2348 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2349 if (m1 != m2) { 2350 coloring->F = 0; 2351 } 2352 } 2353 2354 if (coloring->F) { 2355 w1 = coloring->F; 2356 coloring->F = 0; 2357 } else { 2358 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2359 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2360 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2361 } 2362 2363 /* 2364 Compute all the scale factors and share with other processors 2365 */ 2366 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2367 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2368 for (k=0; k<coloring->ncolors; k++) { 2369 /* 2370 Loop over each column associated with color adding the 2371 perturbation to the vector w3. 2372 */ 2373 for (l=0; l<coloring->ncolumns[k]; l++) { 2374 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2375 dx = xx[col]; 2376 if (dx == 0.0) dx = 1.0; 2377 #if !defined(PETSC_USE_COMPLEX) 2378 if (dx < umin && dx >= 0.0) dx = umin; 2379 else if (dx < 0.0 && dx > -umin) dx = -umin; 2380 #else 2381 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2382 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2383 #endif 2384 dx *= epsilon; 2385 vscale_array[col] = 1.0/dx; 2386 } 2387 } 2388 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2389 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2390 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2391 2392 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2393 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2394 2395 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2396 else vscaleforrow = coloring->columnsforrow; 2397 2398 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2399 /* 2400 Loop over each color 2401 */ 2402 for (k=0; k<coloring->ncolors; k++) { 2403 coloring->currentcolor = k; 2404 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2405 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2406 /* 2407 Loop over each column associated with color adding the 2408 perturbation to the vector w3. 2409 */ 2410 for (l=0; l<coloring->ncolumns[k]; l++) { 2411 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2412 dx = xx[col]; 2413 if (dx == 0.0) dx = 1.0; 2414 #if !defined(PETSC_USE_COMPLEX) 2415 if (dx < umin && dx >= 0.0) dx = umin; 2416 else if (dx < 0.0 && dx > -umin) dx = -umin; 2417 #else 2418 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2419 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2420 #endif 2421 dx *= epsilon; 2422 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2423 w3_array[col] += dx; 2424 } 2425 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2426 2427 /* 2428 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2429 */ 2430 2431 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2432 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2433 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2434 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2435 2436 /* 2437 Loop over rows of vector, putting results into Jacobian matrix 2438 */ 2439 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2440 for (l=0; l<coloring->nrows[k]; l++) { 2441 row = coloring->rows[k][l]; 2442 col = coloring->columnsforrow[k][l]; 2443 y[row] *= vscale_array[vscaleforrow[k][l]]; 2444 srow = row + start; 2445 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2446 } 2447 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2448 } 2449 coloring->currentcolor = k; 2450 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2451 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2452 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2453 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2454 PetscFunctionReturn(0); 2455 } 2456 2457 /* 2458 Computes the number of nonzeros per row needed for preallocation when X and Y 2459 have different nonzero structure. 2460 */ 2461 #undef __FUNCT__ 2462 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2463 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt* nnz) 2464 { 2465 PetscInt i,m=Y->rmap->N; 2466 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2467 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2468 const PetscInt *xi = x->i,*yi = y->i; 2469 2470 PetscFunctionBegin; 2471 /* Set the number of nonzeros in the new matrix */ 2472 for(i=0; i<m; i++) { 2473 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2474 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2475 nnz[i] = 0; 2476 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2477 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2478 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2479 nnz[i]++; 2480 } 2481 for (; k<nzy; k++) nnz[i]++; 2482 } 2483 PetscFunctionReturn(0); 2484 } 2485 2486 #undef __FUNCT__ 2487 #define __FUNCT__ "MatAXPY_SeqAIJ" 2488 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2489 { 2490 PetscErrorCode ierr; 2491 PetscInt i; 2492 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2493 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2494 2495 PetscFunctionBegin; 2496 if (str == SAME_NONZERO_PATTERN) { 2497 PetscScalar alpha = a; 2498 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2499 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2500 if (y->xtoy && y->XtoY != X) { 2501 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2502 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2503 } 2504 if (!y->xtoy) { /* get xtoy */ 2505 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2506 y->XtoY = X; 2507 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2508 } 2509 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2510 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2511 } else { 2512 Mat B; 2513 PetscInt *nnz; 2514 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2515 ierr = MatCreate(((PetscObject)Y)->comm,&B);CHKERRQ(ierr); 2516 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2517 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2518 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2519 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2520 ierr = MatSeqAIJSetPreallocation(B,PETSC_NULL,nnz);CHKERRQ(ierr); 2521 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2522 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2523 ierr = PetscFree(nnz);CHKERRQ(ierr); 2524 } 2525 PetscFunctionReturn(0); 2526 } 2527 2528 #undef __FUNCT__ 2529 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2530 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2531 { 2532 PetscErrorCode ierr; 2533 2534 PetscFunctionBegin; 2535 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2536 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2537 PetscFunctionReturn(0); 2538 } 2539 2540 #undef __FUNCT__ 2541 #define __FUNCT__ "MatConjugate_SeqAIJ" 2542 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2543 { 2544 #if defined(PETSC_USE_COMPLEX) 2545 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2546 PetscInt i,nz; 2547 PetscScalar *a; 2548 2549 PetscFunctionBegin; 2550 nz = aij->nz; 2551 a = aij->a; 2552 for (i=0; i<nz; i++) { 2553 a[i] = PetscConj(a[i]); 2554 } 2555 #else 2556 PetscFunctionBegin; 2557 #endif 2558 PetscFunctionReturn(0); 2559 } 2560 2561 #undef __FUNCT__ 2562 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2563 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2564 { 2565 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2566 PetscErrorCode ierr; 2567 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2568 PetscReal atmp; 2569 PetscScalar *x; 2570 MatScalar *aa; 2571 2572 PetscFunctionBegin; 2573 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2574 aa = a->a; 2575 ai = a->i; 2576 aj = a->j; 2577 2578 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2579 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2580 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2581 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2582 for (i=0; i<m; i++) { 2583 ncols = ai[1] - ai[0]; ai++; 2584 x[i] = 0.0; 2585 for (j=0; j<ncols; j++){ 2586 atmp = PetscAbsScalar(*aa); 2587 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2588 aa++; aj++; 2589 } 2590 } 2591 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2592 PetscFunctionReturn(0); 2593 } 2594 2595 #undef __FUNCT__ 2596 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2597 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2598 { 2599 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2600 PetscErrorCode ierr; 2601 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2602 PetscScalar *x; 2603 MatScalar *aa; 2604 2605 PetscFunctionBegin; 2606 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2607 aa = a->a; 2608 ai = a->i; 2609 aj = a->j; 2610 2611 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2612 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2613 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2614 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2615 for (i=0; i<m; i++) { 2616 ncols = ai[1] - ai[0]; ai++; 2617 if (ncols == A->cmap->n) { /* row is dense */ 2618 x[i] = *aa; if (idx) idx[i] = 0; 2619 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2620 x[i] = 0.0; 2621 if (idx) { 2622 idx[i] = 0; /* in case ncols is zero */ 2623 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2624 if (aj[j] > j) { 2625 idx[i] = j; 2626 break; 2627 } 2628 } 2629 } 2630 } 2631 for (j=0; j<ncols; j++){ 2632 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2633 aa++; aj++; 2634 } 2635 } 2636 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 #undef __FUNCT__ 2641 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2642 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2643 { 2644 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2645 PetscErrorCode ierr; 2646 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2647 PetscReal atmp; 2648 PetscScalar *x; 2649 MatScalar *aa; 2650 2651 PetscFunctionBegin; 2652 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2653 aa = a->a; 2654 ai = a->i; 2655 aj = a->j; 2656 2657 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2658 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2659 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2660 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2661 for (i=0; i<m; i++) { 2662 ncols = ai[1] - ai[0]; ai++; 2663 if (ncols) { 2664 /* Get first nonzero */ 2665 for(j = 0; j < ncols; j++) { 2666 atmp = PetscAbsScalar(aa[j]); 2667 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2668 } 2669 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2670 } else { 2671 x[i] = 0.0; if (idx) idx[i] = 0; 2672 } 2673 for(j = 0; j < ncols; j++) { 2674 atmp = PetscAbsScalar(*aa); 2675 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2676 aa++; aj++; 2677 } 2678 } 2679 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2680 PetscFunctionReturn(0); 2681 } 2682 2683 #undef __FUNCT__ 2684 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2685 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2686 { 2687 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2688 PetscErrorCode ierr; 2689 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2690 PetscScalar *x; 2691 MatScalar *aa; 2692 2693 PetscFunctionBegin; 2694 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2695 aa = a->a; 2696 ai = a->i; 2697 aj = a->j; 2698 2699 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2700 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2701 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2702 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2703 for (i=0; i<m; i++) { 2704 ncols = ai[1] - ai[0]; ai++; 2705 if (ncols == A->cmap->n) { /* row is dense */ 2706 x[i] = *aa; if (idx) idx[i] = 0; 2707 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2708 x[i] = 0.0; 2709 if (idx) { /* find first implicit 0.0 in the row */ 2710 idx[i] = 0; /* in case ncols is zero */ 2711 for (j=0;j<ncols;j++) { 2712 if (aj[j] > j) { 2713 idx[i] = j; 2714 break; 2715 } 2716 } 2717 } 2718 } 2719 for (j=0; j<ncols; j++){ 2720 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2721 aa++; aj++; 2722 } 2723 } 2724 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2725 PetscFunctionReturn(0); 2726 } 2727 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2728 /* -------------------------------------------------------------------*/ 2729 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2730 MatGetRow_SeqAIJ, 2731 MatRestoreRow_SeqAIJ, 2732 MatMult_SeqAIJ, 2733 /* 4*/ MatMultAdd_SeqAIJ, 2734 MatMultTranspose_SeqAIJ, 2735 MatMultTransposeAdd_SeqAIJ, 2736 0, 2737 0, 2738 0, 2739 /*10*/ 0, 2740 MatLUFactor_SeqAIJ, 2741 0, 2742 MatSOR_SeqAIJ, 2743 MatTranspose_SeqAIJ, 2744 /*15*/ MatGetInfo_SeqAIJ, 2745 MatEqual_SeqAIJ, 2746 MatGetDiagonal_SeqAIJ, 2747 MatDiagonalScale_SeqAIJ, 2748 MatNorm_SeqAIJ, 2749 /*20*/ 0, 2750 MatAssemblyEnd_SeqAIJ, 2751 MatSetOption_SeqAIJ, 2752 MatZeroEntries_SeqAIJ, 2753 /*24*/ MatZeroRows_SeqAIJ, 2754 0, 2755 0, 2756 0, 2757 0, 2758 /*29*/ MatSetUpPreallocation_SeqAIJ, 2759 0, 2760 0, 2761 MatGetArray_SeqAIJ, 2762 MatRestoreArray_SeqAIJ, 2763 /*34*/ MatDuplicate_SeqAIJ, 2764 0, 2765 0, 2766 MatILUFactor_SeqAIJ, 2767 0, 2768 /*39*/ MatAXPY_SeqAIJ, 2769 MatGetSubMatrices_SeqAIJ, 2770 MatIncreaseOverlap_SeqAIJ, 2771 MatGetValues_SeqAIJ, 2772 MatCopy_SeqAIJ, 2773 /*44*/ MatGetRowMax_SeqAIJ, 2774 MatScale_SeqAIJ, 2775 0, 2776 MatDiagonalSet_SeqAIJ, 2777 MatZeroRowsColumns_SeqAIJ, 2778 /*49*/ MatSetBlockSize_SeqAIJ, 2779 MatGetRowIJ_SeqAIJ, 2780 MatRestoreRowIJ_SeqAIJ, 2781 MatGetColumnIJ_SeqAIJ, 2782 MatRestoreColumnIJ_SeqAIJ, 2783 /*54*/ MatFDColoringCreate_SeqAIJ, 2784 0, 2785 0, 2786 MatPermute_SeqAIJ, 2787 0, 2788 /*59*/ 0, 2789 MatDestroy_SeqAIJ, 2790 MatView_SeqAIJ, 2791 0, 2792 0, 2793 /*64*/ 0, 2794 0, 2795 0, 2796 0, 2797 0, 2798 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2799 MatGetRowMinAbs_SeqAIJ, 2800 0, 2801 MatSetColoring_SeqAIJ, 2802 #if defined(PETSC_HAVE_ADIC) 2803 MatSetValuesAdic_SeqAIJ, 2804 #else 2805 0, 2806 #endif 2807 /*74*/ MatSetValuesAdifor_SeqAIJ, 2808 MatFDColoringApply_AIJ, 2809 0, 2810 0, 2811 0, 2812 /*79*/ MatFindZeroDiagonals_SeqAIJ, 2813 0, 2814 0, 2815 0, 2816 MatLoad_SeqAIJ, 2817 /*84*/ MatIsSymmetric_SeqAIJ, 2818 MatIsHermitian_SeqAIJ, 2819 0, 2820 0, 2821 0, 2822 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2823 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2824 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2825 MatPtAP_Basic, 2826 MatPtAPSymbolic_SeqAIJ, 2827 /*94*/ MatPtAPNumeric_SeqAIJ, 2828 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2829 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2830 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2831 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2832 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2833 0, 2834 0, 2835 MatConjugate_SeqAIJ, 2836 0, 2837 /*104*/MatSetValuesRow_SeqAIJ, 2838 MatRealPart_SeqAIJ, 2839 MatImaginaryPart_SeqAIJ, 2840 0, 2841 0, 2842 /*109*/MatMatSolve_SeqAIJ, 2843 0, 2844 MatGetRowMin_SeqAIJ, 2845 0, 2846 MatMissingDiagonal_SeqAIJ, 2847 /*114*/0, 2848 0, 2849 0, 2850 0, 2851 0, 2852 /*119*/0, 2853 0, 2854 0, 2855 0, 2856 MatGetMultiProcBlock_SeqAIJ, 2857 /*124*/MatFindNonzeroRows_SeqAIJ, 2858 MatGetColumnNorms_SeqAIJ 2859 }; 2860 2861 EXTERN_C_BEGIN 2862 #undef __FUNCT__ 2863 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2864 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2865 { 2866 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2867 PetscInt i,nz,n; 2868 2869 PetscFunctionBegin; 2870 2871 nz = aij->maxnz; 2872 n = mat->rmap->n; 2873 for (i=0; i<nz; i++) { 2874 aij->j[i] = indices[i]; 2875 } 2876 aij->nz = nz; 2877 for (i=0; i<n; i++) { 2878 aij->ilen[i] = aij->imax[i]; 2879 } 2880 2881 PetscFunctionReturn(0); 2882 } 2883 EXTERN_C_END 2884 2885 #undef __FUNCT__ 2886 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2887 /*@ 2888 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2889 in the matrix. 2890 2891 Input Parameters: 2892 + mat - the SeqAIJ matrix 2893 - indices - the column indices 2894 2895 Level: advanced 2896 2897 Notes: 2898 This can be called if you have precomputed the nonzero structure of the 2899 matrix and want to provide it to the matrix object to improve the performance 2900 of the MatSetValues() operation. 2901 2902 You MUST have set the correct numbers of nonzeros per row in the call to 2903 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2904 2905 MUST be called before any calls to MatSetValues(); 2906 2907 The indices should start with zero, not one. 2908 2909 @*/ 2910 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2911 { 2912 PetscErrorCode ierr; 2913 2914 PetscFunctionBegin; 2915 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2916 PetscValidPointer(indices,2); 2917 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 2918 PetscFunctionReturn(0); 2919 } 2920 2921 /* ----------------------------------------------------------------------------------------*/ 2922 2923 EXTERN_C_BEGIN 2924 #undef __FUNCT__ 2925 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2926 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 2927 { 2928 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2929 PetscErrorCode ierr; 2930 size_t nz = aij->i[mat->rmap->n]; 2931 2932 PetscFunctionBegin; 2933 if (aij->nonew != 1) { 2934 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2935 } 2936 2937 /* allocate space for values if not already there */ 2938 if (!aij->saved_values) { 2939 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2940 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2941 } 2942 2943 /* copy values over */ 2944 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2945 PetscFunctionReturn(0); 2946 } 2947 EXTERN_C_END 2948 2949 #undef __FUNCT__ 2950 #define __FUNCT__ "MatStoreValues" 2951 /*@ 2952 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2953 example, reuse of the linear part of a Jacobian, while recomputing the 2954 nonlinear portion. 2955 2956 Collect on Mat 2957 2958 Input Parameters: 2959 . mat - the matrix (currently only AIJ matrices support this option) 2960 2961 Level: advanced 2962 2963 Common Usage, with SNESSolve(): 2964 $ Create Jacobian matrix 2965 $ Set linear terms into matrix 2966 $ Apply boundary conditions to matrix, at this time matrix must have 2967 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2968 $ boundary conditions again will not change the nonzero structure 2969 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2970 $ ierr = MatStoreValues(mat); 2971 $ Call SNESSetJacobian() with matrix 2972 $ In your Jacobian routine 2973 $ ierr = MatRetrieveValues(mat); 2974 $ Set nonlinear terms in matrix 2975 2976 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2977 $ // build linear portion of Jacobian 2978 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2979 $ ierr = MatStoreValues(mat); 2980 $ loop over nonlinear iterations 2981 $ ierr = MatRetrieveValues(mat); 2982 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2983 $ // call MatAssemblyBegin/End() on matrix 2984 $ Solve linear system with Jacobian 2985 $ endloop 2986 2987 Notes: 2988 Matrix must already be assemblied before calling this routine 2989 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2990 calling this routine. 2991 2992 When this is called multiple times it overwrites the previous set of stored values 2993 and does not allocated additional space. 2994 2995 .seealso: MatRetrieveValues() 2996 2997 @*/ 2998 PetscErrorCode MatStoreValues(Mat mat) 2999 { 3000 PetscErrorCode ierr; 3001 3002 PetscFunctionBegin; 3003 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3004 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3005 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3006 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3007 PetscFunctionReturn(0); 3008 } 3009 3010 EXTERN_C_BEGIN 3011 #undef __FUNCT__ 3012 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3013 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3014 { 3015 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3016 PetscErrorCode ierr; 3017 PetscInt nz = aij->i[mat->rmap->n]; 3018 3019 PetscFunctionBegin; 3020 if (aij->nonew != 1) { 3021 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3022 } 3023 if (!aij->saved_values) { 3024 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3025 } 3026 /* copy values over */ 3027 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3028 PetscFunctionReturn(0); 3029 } 3030 EXTERN_C_END 3031 3032 #undef __FUNCT__ 3033 #define __FUNCT__ "MatRetrieveValues" 3034 /*@ 3035 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3036 example, reuse of the linear part of a Jacobian, while recomputing the 3037 nonlinear portion. 3038 3039 Collect on Mat 3040 3041 Input Parameters: 3042 . mat - the matrix (currently on AIJ matrices support this option) 3043 3044 Level: advanced 3045 3046 .seealso: MatStoreValues() 3047 3048 @*/ 3049 PetscErrorCode MatRetrieveValues(Mat mat) 3050 { 3051 PetscErrorCode ierr; 3052 3053 PetscFunctionBegin; 3054 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3055 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3056 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3057 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3058 PetscFunctionReturn(0); 3059 } 3060 3061 3062 /* --------------------------------------------------------------------------------*/ 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "MatCreateSeqAIJ" 3065 /*@C 3066 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3067 (the default parallel PETSc format). For good matrix assembly performance 3068 the user should preallocate the matrix storage by setting the parameter nz 3069 (or the array nnz). By setting these parameters accurately, performance 3070 during matrix assembly can be increased by more than a factor of 50. 3071 3072 Collective on MPI_Comm 3073 3074 Input Parameters: 3075 + comm - MPI communicator, set to PETSC_COMM_SELF 3076 . m - number of rows 3077 . n - number of columns 3078 . nz - number of nonzeros per row (same for all rows) 3079 - nnz - array containing the number of nonzeros in the various rows 3080 (possibly different for each row) or PETSC_NULL 3081 3082 Output Parameter: 3083 . A - the matrix 3084 3085 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3086 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3087 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3088 3089 Notes: 3090 If nnz is given then nz is ignored 3091 3092 The AIJ format (also called the Yale sparse matrix format or 3093 compressed row storage), is fully compatible with standard Fortran 77 3094 storage. That is, the stored row and column indices can begin at 3095 either one (as in Fortran) or zero. See the users' manual for details. 3096 3097 Specify the preallocated storage with either nz or nnz (not both). 3098 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3099 allocation. For large problems you MUST preallocate memory or you 3100 will get TERRIBLE performance, see the users' manual chapter on matrices. 3101 3102 By default, this format uses inodes (identical nodes) when possible, to 3103 improve numerical efficiency of matrix-vector products and solves. We 3104 search for consecutive rows with the same nonzero structure, thereby 3105 reusing matrix information to achieve increased efficiency. 3106 3107 Options Database Keys: 3108 + -mat_no_inode - Do not use inodes 3109 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3110 3111 Level: intermediate 3112 3113 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3114 3115 @*/ 3116 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3117 { 3118 PetscErrorCode ierr; 3119 3120 PetscFunctionBegin; 3121 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3122 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3123 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3124 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3125 PetscFunctionReturn(0); 3126 } 3127 3128 #undef __FUNCT__ 3129 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3130 /*@C 3131 MatSeqAIJSetPreallocation - For good matrix assembly performance 3132 the user should preallocate the matrix storage by setting the parameter nz 3133 (or the array nnz). By setting these parameters accurately, performance 3134 during matrix assembly can be increased by more than a factor of 50. 3135 3136 Collective on MPI_Comm 3137 3138 Input Parameters: 3139 + B - The matrix-free 3140 . nz - number of nonzeros per row (same for all rows) 3141 - nnz - array containing the number of nonzeros in the various rows 3142 (possibly different for each row) or PETSC_NULL 3143 3144 Notes: 3145 If nnz is given then nz is ignored 3146 3147 The AIJ format (also called the Yale sparse matrix format or 3148 compressed row storage), is fully compatible with standard Fortran 77 3149 storage. That is, the stored row and column indices can begin at 3150 either one (as in Fortran) or zero. See the users' manual for details. 3151 3152 Specify the preallocated storage with either nz or nnz (not both). 3153 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3154 allocation. For large problems you MUST preallocate memory or you 3155 will get TERRIBLE performance, see the users' manual chapter on matrices. 3156 3157 You can call MatGetInfo() to get information on how effective the preallocation was; 3158 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3159 You can also run with the option -info and look for messages with the string 3160 malloc in them to see if additional memory allocation was needed. 3161 3162 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3163 entries or columns indices 3164 3165 By default, this format uses inodes (identical nodes) when possible, to 3166 improve numerical efficiency of matrix-vector products and solves. We 3167 search for consecutive rows with the same nonzero structure, thereby 3168 reusing matrix information to achieve increased efficiency. 3169 3170 Options Database Keys: 3171 + -mat_no_inode - Do not use inodes 3172 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3173 - -mat_aij_oneindex - Internally use indexing starting at 1 3174 rather than 0. Note that when calling MatSetValues(), 3175 the user still MUST index entries starting at 0! 3176 3177 Level: intermediate 3178 3179 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3180 3181 @*/ 3182 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3183 { 3184 PetscErrorCode ierr; 3185 3186 PetscFunctionBegin; 3187 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3188 PetscFunctionReturn(0); 3189 } 3190 3191 EXTERN_C_BEGIN 3192 #undef __FUNCT__ 3193 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3194 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3195 { 3196 Mat_SeqAIJ *b; 3197 PetscBool skipallocation = PETSC_FALSE; 3198 PetscErrorCode ierr; 3199 PetscInt i; 3200 3201 PetscFunctionBegin; 3202 3203 if (nz == MAT_SKIP_ALLOCATION) { 3204 skipallocation = PETSC_TRUE; 3205 nz = 0; 3206 } 3207 3208 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3209 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3210 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3211 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3212 3213 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3214 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3215 if (nnz) { 3216 for (i=0; i<B->rmap->n; i++) { 3217 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]); 3218 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); 3219 } 3220 } 3221 3222 B->preallocated = PETSC_TRUE; 3223 b = (Mat_SeqAIJ*)B->data; 3224 3225 if (!skipallocation) { 3226 if (!b->imax) { 3227 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3228 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3229 } 3230 if (!nnz) { 3231 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3232 else if (nz < 0) nz = 1; 3233 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3234 nz = nz*B->rmap->n; 3235 } else { 3236 nz = 0; 3237 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3238 } 3239 /* b->ilen will count nonzeros in each row so far. */ 3240 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3241 3242 /* allocate the matrix space */ 3243 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3244 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3245 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3246 b->i[0] = 0; 3247 for (i=1; i<B->rmap->n+1; i++) { 3248 b->i[i] = b->i[i-1] + b->imax[i-1]; 3249 } 3250 b->singlemalloc = PETSC_TRUE; 3251 b->free_a = PETSC_TRUE; 3252 b->free_ij = PETSC_TRUE; 3253 } else { 3254 b->free_a = PETSC_FALSE; 3255 b->free_ij = PETSC_FALSE; 3256 } 3257 3258 b->nz = 0; 3259 b->maxnz = nz; 3260 B->info.nz_unneeded = (double)b->maxnz; 3261 PetscFunctionReturn(0); 3262 } 3263 EXTERN_C_END 3264 3265 #undef __FUNCT__ 3266 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3267 /*@ 3268 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3269 3270 Input Parameters: 3271 + B - the matrix 3272 . i - the indices into j for the start of each row (starts with zero) 3273 . j - the column indices for each row (starts with zero) these must be sorted for each row 3274 - v - optional values in the matrix 3275 3276 Level: developer 3277 3278 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3279 3280 .keywords: matrix, aij, compressed row, sparse, sequential 3281 3282 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3283 @*/ 3284 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3285 { 3286 PetscErrorCode ierr; 3287 3288 PetscFunctionBegin; 3289 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3290 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3291 PetscFunctionReturn(0); 3292 } 3293 3294 EXTERN_C_BEGIN 3295 #undef __FUNCT__ 3296 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3297 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3298 { 3299 PetscInt i; 3300 PetscInt m,n; 3301 PetscInt nz; 3302 PetscInt *nnz, nz_max = 0; 3303 PetscScalar *values; 3304 PetscErrorCode ierr; 3305 3306 PetscFunctionBegin; 3307 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3308 3309 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3310 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3311 for(i = 0; i < m; i++) { 3312 nz = Ii[i+1]- Ii[i]; 3313 nz_max = PetscMax(nz_max, nz); 3314 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3315 nnz[i] = nz; 3316 } 3317 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3318 ierr = PetscFree(nnz);CHKERRQ(ierr); 3319 3320 if (v) { 3321 values = (PetscScalar*) v; 3322 } else { 3323 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3324 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3325 } 3326 3327 for(i = 0; i < m; i++) { 3328 nz = Ii[i+1] - Ii[i]; 3329 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3330 } 3331 3332 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3333 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3334 3335 if (!v) { 3336 ierr = PetscFree(values);CHKERRQ(ierr); 3337 } 3338 PetscFunctionReturn(0); 3339 } 3340 EXTERN_C_END 3341 3342 #include <../src/mat/impls/dense/seq/dense.h> 3343 #include <private/petscaxpy.h> 3344 3345 #undef __FUNCT__ 3346 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3347 /* 3348 Computes (B'*A')' since computing B*A directly is untenable 3349 3350 n p p 3351 ( ) ( ) ( ) 3352 m ( A ) * n ( B ) = m ( C ) 3353 ( ) ( ) ( ) 3354 3355 */ 3356 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3357 { 3358 PetscErrorCode ierr; 3359 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3360 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3361 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3362 PetscInt i,n,m,q,p; 3363 const PetscInt *ii,*idx; 3364 const PetscScalar *b,*a,*a_q; 3365 PetscScalar *c,*c_q; 3366 3367 PetscFunctionBegin; 3368 m = A->rmap->n; 3369 n = A->cmap->n; 3370 p = B->cmap->n; 3371 a = sub_a->v; 3372 b = sub_b->a; 3373 c = sub_c->v; 3374 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3375 3376 ii = sub_b->i; 3377 idx = sub_b->j; 3378 for (i=0; i<n; i++) { 3379 q = ii[i+1] - ii[i]; 3380 while (q-->0) { 3381 c_q = c + m*(*idx); 3382 a_q = a + m*i; 3383 PetscAXPY(c_q,*b,a_q,m); 3384 idx++; 3385 b++; 3386 } 3387 } 3388 PetscFunctionReturn(0); 3389 } 3390 3391 #undef __FUNCT__ 3392 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3393 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3394 { 3395 PetscErrorCode ierr; 3396 PetscInt m=A->rmap->n,n=B->cmap->n; 3397 Mat Cmat; 3398 3399 PetscFunctionBegin; 3400 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); 3401 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3402 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3403 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3404 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3405 Cmat->assembled = PETSC_TRUE; 3406 *C = Cmat; 3407 PetscFunctionReturn(0); 3408 } 3409 3410 /* ----------------------------------------------------------------*/ 3411 #undef __FUNCT__ 3412 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3413 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3414 { 3415 PetscErrorCode ierr; 3416 3417 PetscFunctionBegin; 3418 if (scall == MAT_INITIAL_MATRIX){ 3419 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3420 } 3421 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3422 PetscFunctionReturn(0); 3423 } 3424 3425 3426 /*MC 3427 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3428 based on compressed sparse row format. 3429 3430 Options Database Keys: 3431 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3432 3433 Level: beginner 3434 3435 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3436 M*/ 3437 3438 EXTERN_C_BEGIN 3439 #if defined(PETSC_HAVE_PASTIX) 3440 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3441 #endif 3442 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3443 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3444 #endif 3445 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3446 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3447 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3448 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3449 #if defined(PETSC_HAVE_MUMPS) 3450 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3451 #endif 3452 #if defined(PETSC_HAVE_SUPERLU) 3453 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3454 #endif 3455 #if defined(PETSC_HAVE_SUPERLU_DIST) 3456 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3457 #endif 3458 #if defined(PETSC_HAVE_SPOOLES) 3459 extern PetscErrorCode MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3460 #endif 3461 #if defined(PETSC_HAVE_UMFPACK) 3462 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3463 #endif 3464 #if defined(PETSC_HAVE_CHOLMOD) 3465 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3466 #endif 3467 #if defined(PETSC_HAVE_LUSOL) 3468 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3469 #endif 3470 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3471 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3472 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3473 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3474 #endif 3475 EXTERN_C_END 3476 3477 EXTERN_C_BEGIN 3478 #undef __FUNCT__ 3479 #define __FUNCT__ "MatCreate_SeqAIJ" 3480 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3481 { 3482 Mat_SeqAIJ *b; 3483 PetscErrorCode ierr; 3484 PetscMPIInt size; 3485 3486 PetscFunctionBegin; 3487 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3488 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3489 3490 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3491 B->data = (void*)b; 3492 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3493 b->row = 0; 3494 b->col = 0; 3495 b->icol = 0; 3496 b->reallocs = 0; 3497 b->ignorezeroentries = PETSC_FALSE; 3498 b->roworiented = PETSC_TRUE; 3499 b->nonew = 0; 3500 b->diag = 0; 3501 b->solve_work = 0; 3502 B->spptr = 0; 3503 b->saved_values = 0; 3504 b->idiag = 0; 3505 b->mdiag = 0; 3506 b->ssor_work = 0; 3507 b->omega = 1.0; 3508 b->fshift = 0.0; 3509 b->idiagvalid = PETSC_FALSE; 3510 b->keepnonzeropattern = PETSC_FALSE; 3511 b->xtoy = 0; 3512 b->XtoY = 0; 3513 B->same_nonzero = PETSC_FALSE; 3514 3515 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3516 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3517 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3518 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3519 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3520 #endif 3521 #if defined(PETSC_HAVE_PASTIX) 3522 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3523 #endif 3524 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3525 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3526 #endif 3527 #if defined(PETSC_HAVE_SUPERLU) 3528 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3529 #endif 3530 #if defined(PETSC_HAVE_SUPERLU_DIST) 3531 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3532 #endif 3533 #if defined(PETSC_HAVE_SPOOLES) 3534 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C","MatGetFactor_seqaij_spooles",MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3535 #endif 3536 #if defined(PETSC_HAVE_MUMPS) 3537 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3538 #endif 3539 #if defined(PETSC_HAVE_UMFPACK) 3540 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3541 #endif 3542 #if defined(PETSC_HAVE_CHOLMOD) 3543 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3544 #endif 3545 #if defined(PETSC_HAVE_LUSOL) 3546 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_seqaij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3547 #endif 3548 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3549 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3550 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3551 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3552 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3553 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3554 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3555 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3556 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3557 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3558 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3559 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3560 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3561 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3562 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3563 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3564 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3565 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3566 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3567 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3568 PetscFunctionReturn(0); 3569 } 3570 EXTERN_C_END 3571 3572 #undef __FUNCT__ 3573 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3574 /* 3575 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3576 */ 3577 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3578 { 3579 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3580 PetscErrorCode ierr; 3581 PetscInt i,m = A->rmap->n; 3582 3583 PetscFunctionBegin; 3584 c = (Mat_SeqAIJ*)C->data; 3585 3586 C->factortype = A->factortype; 3587 c->row = 0; 3588 c->col = 0; 3589 c->icol = 0; 3590 c->reallocs = 0; 3591 3592 C->assembled = PETSC_TRUE; 3593 3594 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3595 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3596 3597 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3598 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3599 for (i=0; i<m; i++) { 3600 c->imax[i] = a->imax[i]; 3601 c->ilen[i] = a->ilen[i]; 3602 } 3603 3604 /* allocate the matrix space */ 3605 if (mallocmatspace){ 3606 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3607 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3608 c->singlemalloc = PETSC_TRUE; 3609 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3610 if (m > 0) { 3611 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3612 if (cpvalues == MAT_COPY_VALUES) { 3613 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3614 } else { 3615 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3616 } 3617 } 3618 } 3619 3620 c->ignorezeroentries = a->ignorezeroentries; 3621 c->roworiented = a->roworiented; 3622 c->nonew = a->nonew; 3623 if (a->diag) { 3624 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3625 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3626 for (i=0; i<m; i++) { 3627 c->diag[i] = a->diag[i]; 3628 } 3629 } else c->diag = 0; 3630 c->solve_work = 0; 3631 c->saved_values = 0; 3632 c->idiag = 0; 3633 c->ssor_work = 0; 3634 c->keepnonzeropattern = a->keepnonzeropattern; 3635 c->free_a = PETSC_TRUE; 3636 c->free_ij = PETSC_TRUE; 3637 c->xtoy = 0; 3638 c->XtoY = 0; 3639 3640 c->nz = a->nz; 3641 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3642 C->preallocated = PETSC_TRUE; 3643 3644 c->compressedrow.use = a->compressedrow.use; 3645 c->compressedrow.nrows = a->compressedrow.nrows; 3646 c->compressedrow.check = a->compressedrow.check; 3647 if (a->compressedrow.use){ 3648 i = a->compressedrow.nrows; 3649 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3650 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3651 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3652 } else { 3653 c->compressedrow.use = PETSC_FALSE; 3654 c->compressedrow.i = PETSC_NULL; 3655 c->compressedrow.rindex = PETSC_NULL; 3656 } 3657 C->same_nonzero = A->same_nonzero; 3658 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3659 3660 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3661 PetscFunctionReturn(0); 3662 } 3663 3664 #undef __FUNCT__ 3665 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3666 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3667 { 3668 PetscErrorCode ierr; 3669 3670 PetscFunctionBegin; 3671 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3672 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3673 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3674 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3675 PetscFunctionReturn(0); 3676 } 3677 3678 #undef __FUNCT__ 3679 #define __FUNCT__ "MatLoad_SeqAIJ" 3680 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3681 { 3682 Mat_SeqAIJ *a; 3683 PetscErrorCode ierr; 3684 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3685 int fd; 3686 PetscMPIInt size; 3687 MPI_Comm comm; 3688 3689 PetscFunctionBegin; 3690 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3691 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3692 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3693 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3694 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3695 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3696 M = header[1]; N = header[2]; nz = header[3]; 3697 3698 if (nz < 0) { 3699 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3700 } 3701 3702 /* read in row lengths */ 3703 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3704 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3705 3706 /* check if sum of rowlengths is same as nz */ 3707 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3708 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); 3709 3710 /* set global size if not set already*/ 3711 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3712 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3713 } else { 3714 /* if sizes and type are already set, check if the vector global sizes are correct */ 3715 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3716 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); 3717 } 3718 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3719 a = (Mat_SeqAIJ*)newMat->data; 3720 3721 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3722 3723 /* read in nonzero values */ 3724 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3725 3726 /* set matrix "i" values */ 3727 a->i[0] = 0; 3728 for (i=1; i<= M; i++) { 3729 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3730 a->ilen[i-1] = rowlengths[i-1]; 3731 } 3732 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3733 3734 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3735 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3736 PetscFunctionReturn(0); 3737 } 3738 3739 #undef __FUNCT__ 3740 #define __FUNCT__ "MatEqual_SeqAIJ" 3741 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 3742 { 3743 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3744 PetscErrorCode ierr; 3745 #if defined(PETSC_USE_COMPLEX) 3746 PetscInt k; 3747 #endif 3748 3749 PetscFunctionBegin; 3750 /* If the matrix dimensions are not equal,or no of nonzeros */ 3751 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3752 *flg = PETSC_FALSE; 3753 PetscFunctionReturn(0); 3754 } 3755 3756 /* if the a->i are the same */ 3757 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3758 if (!*flg) PetscFunctionReturn(0); 3759 3760 /* if a->j are the same */ 3761 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3762 if (!*flg) PetscFunctionReturn(0); 3763 3764 /* if a->a are the same */ 3765 #if defined(PETSC_USE_COMPLEX) 3766 for (k=0; k<a->nz; k++){ 3767 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3768 *flg = PETSC_FALSE; 3769 PetscFunctionReturn(0); 3770 } 3771 } 3772 #else 3773 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3774 #endif 3775 PetscFunctionReturn(0); 3776 } 3777 3778 #undef __FUNCT__ 3779 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3780 /*@ 3781 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3782 provided by the user. 3783 3784 Collective on MPI_Comm 3785 3786 Input Parameters: 3787 + comm - must be an MPI communicator of size 1 3788 . m - number of rows 3789 . n - number of columns 3790 . i - row indices 3791 . j - column indices 3792 - a - matrix values 3793 3794 Output Parameter: 3795 . mat - the matrix 3796 3797 Level: intermediate 3798 3799 Notes: 3800 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3801 once the matrix is destroyed and not before 3802 3803 You cannot set new nonzero locations into this matrix, that will generate an error. 3804 3805 The i and j indices are 0 based 3806 3807 The format which is used for the sparse matrix input, is equivalent to a 3808 row-major ordering.. i.e for the following matrix, the input data expected is 3809 as shown: 3810 3811 1 0 0 3812 2 0 3 3813 4 5 6 3814 3815 i = {0,1,3,6} [size = nrow+1 = 3+1] 3816 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3817 v = {1,2,3,4,5,6} [size = nz = 6] 3818 3819 3820 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3821 3822 @*/ 3823 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3824 { 3825 PetscErrorCode ierr; 3826 PetscInt ii; 3827 Mat_SeqAIJ *aij; 3828 #if defined(PETSC_USE_DEBUG) 3829 PetscInt jj; 3830 #endif 3831 3832 PetscFunctionBegin; 3833 if (i[0]) { 3834 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3835 } 3836 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3837 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3838 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3839 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3840 aij = (Mat_SeqAIJ*)(*mat)->data; 3841 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3842 3843 aij->i = i; 3844 aij->j = j; 3845 aij->a = a; 3846 aij->singlemalloc = PETSC_FALSE; 3847 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3848 aij->free_a = PETSC_FALSE; 3849 aij->free_ij = PETSC_FALSE; 3850 3851 for (ii=0; ii<m; ii++) { 3852 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3853 #if defined(PETSC_USE_DEBUG) 3854 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]); 3855 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3856 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); 3857 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); 3858 } 3859 #endif 3860 } 3861 #if defined(PETSC_USE_DEBUG) 3862 for (ii=0; ii<aij->i[m]; ii++) { 3863 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3864 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]); 3865 } 3866 #endif 3867 3868 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3869 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3870 PetscFunctionReturn(0); 3871 } 3872 3873 #undef __FUNCT__ 3874 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3875 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3876 { 3877 PetscErrorCode ierr; 3878 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3879 3880 PetscFunctionBegin; 3881 if (coloring->ctype == IS_COLORING_GLOBAL) { 3882 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3883 a->coloring = coloring; 3884 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3885 PetscInt i,*larray; 3886 ISColoring ocoloring; 3887 ISColoringValue *colors; 3888 3889 /* set coloring for diagonal portion */ 3890 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3891 for (i=0; i<A->cmap->n; i++) { 3892 larray[i] = i; 3893 } 3894 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3895 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3896 for (i=0; i<A->cmap->n; i++) { 3897 colors[i] = coloring->colors[larray[i]]; 3898 } 3899 ierr = PetscFree(larray);CHKERRQ(ierr); 3900 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3901 a->coloring = ocoloring; 3902 } 3903 PetscFunctionReturn(0); 3904 } 3905 3906 #if defined(PETSC_HAVE_ADIC) 3907 EXTERN_C_BEGIN 3908 #include <adic/ad_utils.h> 3909 EXTERN_C_END 3910 3911 #undef __FUNCT__ 3912 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3913 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3914 { 3915 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3916 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3917 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3918 ISColoringValue *color; 3919 3920 PetscFunctionBegin; 3921 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3922 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3923 color = a->coloring->colors; 3924 /* loop over rows */ 3925 for (i=0; i<m; i++) { 3926 nz = ii[i+1] - ii[i]; 3927 /* loop over columns putting computed value into matrix */ 3928 for (j=0; j<nz; j++) { 3929 *v++ = values[color[*jj++]]; 3930 } 3931 values += nlen; /* jump to next row of derivatives */ 3932 } 3933 PetscFunctionReturn(0); 3934 } 3935 #endif 3936 3937 #undef __FUNCT__ 3938 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3939 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3940 { 3941 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3942 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3943 MatScalar *v = a->a; 3944 PetscScalar *values = (PetscScalar *)advalues; 3945 ISColoringValue *color; 3946 3947 PetscFunctionBegin; 3948 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3949 color = a->coloring->colors; 3950 /* loop over rows */ 3951 for (i=0; i<m; i++) { 3952 nz = ii[i+1] - ii[i]; 3953 /* loop over columns putting computed value into matrix */ 3954 for (j=0; j<nz; j++) { 3955 *v++ = values[color[*jj++]]; 3956 } 3957 values += nl; /* jump to next row of derivatives */ 3958 } 3959 PetscFunctionReturn(0); 3960 } 3961 3962 /* 3963 Special version for direct calls from Fortran 3964 */ 3965 #include <private/fortranimpl.h> 3966 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3967 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3968 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3969 #define matsetvaluesseqaij_ matsetvaluesseqaij 3970 #endif 3971 3972 /* Change these macros so can be used in void function */ 3973 #undef CHKERRQ 3974 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3975 #undef SETERRQ2 3976 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3977 3978 EXTERN_C_BEGIN 3979 #undef __FUNCT__ 3980 #define __FUNCT__ "matsetvaluesseqaij_" 3981 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3982 { 3983 Mat A = *AA; 3984 PetscInt m = *mm, n = *nn; 3985 InsertMode is = *isis; 3986 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3987 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3988 PetscInt *imax,*ai,*ailen; 3989 PetscErrorCode ierr; 3990 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3991 MatScalar *ap,value,*aa; 3992 PetscBool ignorezeroentries = a->ignorezeroentries; 3993 PetscBool roworiented = a->roworiented; 3994 3995 PetscFunctionBegin; 3996 ierr = MatPreallocated(A);CHKERRQ(ierr); 3997 imax = a->imax; 3998 ai = a->i; 3999 ailen = a->ilen; 4000 aj = a->j; 4001 aa = a->a; 4002 4003 for (k=0; k<m; k++) { /* loop over added rows */ 4004 row = im[k]; 4005 if (row < 0) continue; 4006 #if defined(PETSC_USE_DEBUG) 4007 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4008 #endif 4009 rp = aj + ai[row]; ap = aa + ai[row]; 4010 rmax = imax[row]; nrow = ailen[row]; 4011 low = 0; 4012 high = nrow; 4013 for (l=0; l<n; l++) { /* loop over added columns */ 4014 if (in[l] < 0) continue; 4015 #if defined(PETSC_USE_DEBUG) 4016 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4017 #endif 4018 col = in[l]; 4019 if (roworiented) { 4020 value = v[l + k*n]; 4021 } else { 4022 value = v[k + l*m]; 4023 } 4024 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4025 4026 if (col <= lastcol) low = 0; else high = nrow; 4027 lastcol = col; 4028 while (high-low > 5) { 4029 t = (low+high)/2; 4030 if (rp[t] > col) high = t; 4031 else low = t; 4032 } 4033 for (i=low; i<high; i++) { 4034 if (rp[i] > col) break; 4035 if (rp[i] == col) { 4036 if (is == ADD_VALUES) ap[i] += value; 4037 else ap[i] = value; 4038 goto noinsert; 4039 } 4040 } 4041 if (value == 0.0 && ignorezeroentries) goto noinsert; 4042 if (nonew == 1) goto noinsert; 4043 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4044 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4045 N = nrow++ - 1; a->nz++; high++; 4046 /* shift up all the later entries in this row */ 4047 for (ii=N; ii>=i; ii--) { 4048 rp[ii+1] = rp[ii]; 4049 ap[ii+1] = ap[ii]; 4050 } 4051 rp[i] = col; 4052 ap[i] = value; 4053 noinsert:; 4054 low = i + 1; 4055 } 4056 ailen[row] = nrow; 4057 } 4058 A->same_nonzero = PETSC_FALSE; 4059 PetscFunctionReturnVoid(); 4060 } 4061 EXTERN_C_END 4062 4063