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