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