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