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