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,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,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,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,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_LUSOL) 1836 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1837 #endif 1838 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1839 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1840 #endif 1841 PetscFunctionReturn(0); 1842 } 1843 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1844 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1845 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatFactorInfo*,IS,IS,Mat*); 1846 #undef __FUNCT__ 1847 #define __FUNCT__ "MatCopy_SeqAIJ" 1848 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1849 { 1850 int ierr; 1851 1852 PetscFunctionBegin; 1853 /* If the two matrices have the same copy implementation, use fast copy. */ 1854 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1855 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1856 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1857 1858 if (a->i[A->m] != b->i[B->m]) { 1859 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1860 } 1861 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1862 } else { 1863 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1864 } 1865 PetscFunctionReturn(0); 1866 } 1867 1868 #undef __FUNCT__ 1869 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1870 int MatSetUpPreallocation_SeqAIJ(Mat A) 1871 { 1872 int ierr; 1873 1874 PetscFunctionBegin; 1875 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1876 PetscFunctionReturn(0); 1877 } 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "MatGetArray_SeqAIJ" 1881 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1882 { 1883 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1884 PetscFunctionBegin; 1885 *array = a->a; 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1891 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1892 { 1893 PetscFunctionBegin; 1894 PetscFunctionReturn(0); 1895 } 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1899 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1900 { 1901 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1902 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1903 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1904 PetscScalar *vscale_array; 1905 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1906 Vec w1,w2,w3; 1907 void *fctx = coloring->fctx; 1908 PetscTruth flg; 1909 1910 PetscFunctionBegin; 1911 if (!coloring->w1) { 1912 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1913 PetscLogObjectParent(coloring,coloring->w1); 1914 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1915 PetscLogObjectParent(coloring,coloring->w2); 1916 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1917 PetscLogObjectParent(coloring,coloring->w3); 1918 } 1919 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1920 1921 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1922 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1923 if (flg) { 1924 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1925 } else { 1926 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1927 } 1928 1929 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1930 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1931 1932 /* 1933 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1934 coloring->F for the coarser grids from the finest 1935 */ 1936 if (coloring->F) { 1937 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1938 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1939 if (m1 != m2) { 1940 coloring->F = 0; 1941 } 1942 } 1943 1944 if (coloring->F) { 1945 w1 = coloring->F; 1946 coloring->F = 0; 1947 } else { 1948 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1949 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1950 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1951 } 1952 1953 /* 1954 Compute all the scale factors and share with other processors 1955 */ 1956 ierr = VecGetArrayFast(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1957 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1958 for (k=0; k<coloring->ncolors; k++) { 1959 /* 1960 Loop over each column associated with color adding the 1961 perturbation to the vector w3. 1962 */ 1963 for (l=0; l<coloring->ncolumns[k]; l++) { 1964 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1965 dx = xx[col]; 1966 if (dx == 0.0) dx = 1.0; 1967 #if !defined(PETSC_USE_COMPLEX) 1968 if (dx < umin && dx >= 0.0) dx = umin; 1969 else if (dx < 0.0 && dx > -umin) dx = -umin; 1970 #else 1971 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1972 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1973 #endif 1974 dx *= epsilon; 1975 vscale_array[col] = 1.0/dx; 1976 } 1977 } 1978 vscale_array = vscale_array + start;ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1979 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1980 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1981 1982 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1983 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1984 1985 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1986 else vscaleforrow = coloring->columnsforrow; 1987 1988 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1989 /* 1990 Loop over each color 1991 */ 1992 for (k=0; k<coloring->ncolors; k++) { 1993 coloring->currentcolor = k; 1994 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1995 ierr = VecGetArrayFast(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1996 /* 1997 Loop over each column associated with color adding the 1998 perturbation to the vector w3. 1999 */ 2000 for (l=0; l<coloring->ncolumns[k]; l++) { 2001 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2002 dx = xx[col]; 2003 if (dx == 0.0) dx = 1.0; 2004 #if !defined(PETSC_USE_COMPLEX) 2005 if (dx < umin && dx >= 0.0) dx = umin; 2006 else if (dx < 0.0 && dx > -umin) dx = -umin; 2007 #else 2008 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2009 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2010 #endif 2011 dx *= epsilon; 2012 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2013 w3_array[col] += dx; 2014 } 2015 w3_array = w3_array + start; ierr = VecRestoreArrayFast(w3,&w3_array);CHKERRQ(ierr); 2016 2017 /* 2018 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2019 */ 2020 2021 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2022 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2023 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2024 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2025 2026 /* 2027 Loop over rows of vector, putting results into Jacobian matrix 2028 */ 2029 ierr = VecGetArrayFast(w2,&y);CHKERRQ(ierr); 2030 for (l=0; l<coloring->nrows[k]; l++) { 2031 row = coloring->rows[k][l]; 2032 col = coloring->columnsforrow[k][l]; 2033 y[row] *= vscale_array[vscaleforrow[k][l]]; 2034 srow = row + start; 2035 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2036 } 2037 ierr = VecRestoreArrayFast(w2,&y);CHKERRQ(ierr); 2038 } 2039 coloring->currentcolor = k; 2040 ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2041 xx = xx + start; ierr = VecRestoreArrayFast(x1,&xx);CHKERRQ(ierr); 2042 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2043 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #include "petscblaslapack.h" 2048 #undef __FUNCT__ 2049 #define __FUNCT__ "MatAXPY_SeqAIJ" 2050 int MatAXPY_SeqAIJ(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 2051 { 2052 int ierr,one=1,i; 2053 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2054 2055 PetscFunctionBegin; 2056 if (str == SAME_NONZERO_PATTERN) { 2057 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2058 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2059 if (y->xtoy && y->XtoY != X) { 2060 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2061 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2062 } 2063 if (!y->xtoy) { /* get xtoy */ 2064 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2065 y->XtoY = X; 2066 } 2067 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2068 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2069 } else { 2070 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2071 } 2072 PetscFunctionReturn(0); 2073 } 2074 2075 /* -------------------------------------------------------------------*/ 2076 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2077 MatGetRow_SeqAIJ, 2078 MatRestoreRow_SeqAIJ, 2079 MatMult_SeqAIJ, 2080 MatMultAdd_SeqAIJ, 2081 MatMultTranspose_SeqAIJ, 2082 MatMultTransposeAdd_SeqAIJ, 2083 MatSolve_SeqAIJ, 2084 MatSolveAdd_SeqAIJ, 2085 MatSolveTranspose_SeqAIJ, 2086 MatSolveTransposeAdd_SeqAIJ, 2087 MatLUFactor_SeqAIJ, 2088 0, 2089 MatRelax_SeqAIJ, 2090 MatTranspose_SeqAIJ, 2091 MatGetInfo_SeqAIJ, 2092 MatEqual_SeqAIJ, 2093 MatGetDiagonal_SeqAIJ, 2094 MatDiagonalScale_SeqAIJ, 2095 MatNorm_SeqAIJ, 2096 0, 2097 MatAssemblyEnd_SeqAIJ, 2098 MatCompress_SeqAIJ, 2099 MatSetOption_SeqAIJ, 2100 MatZeroEntries_SeqAIJ, 2101 MatZeroRows_SeqAIJ, 2102 MatLUFactorSymbolic_SeqAIJ, 2103 MatLUFactorNumeric_SeqAIJ, 2104 MatCholeskyFactorSymbolic_SeqAIJ, 2105 MatCholeskyFactorNumeric_SeqAIJ, 2106 MatSetUpPreallocation_SeqAIJ, 2107 MatILUFactorSymbolic_SeqAIJ, 2108 MatICCFactorSymbolic_SeqAIJ, 2109 MatGetArray_SeqAIJ, 2110 MatRestoreArray_SeqAIJ, 2111 MatDuplicate_SeqAIJ, 2112 0, 2113 0, 2114 MatILUFactor_SeqAIJ, 2115 0, 2116 MatAXPY_SeqAIJ, 2117 MatGetSubMatrices_SeqAIJ, 2118 MatIncreaseOverlap_SeqAIJ, 2119 MatGetValues_SeqAIJ, 2120 MatCopy_SeqAIJ, 2121 MatPrintHelp_SeqAIJ, 2122 MatScale_SeqAIJ, 2123 0, 2124 0, 2125 MatILUDTFactor_SeqAIJ, 2126 MatGetBlockSize_SeqAIJ, 2127 MatGetRowIJ_SeqAIJ, 2128 MatRestoreRowIJ_SeqAIJ, 2129 MatGetColumnIJ_SeqAIJ, 2130 MatRestoreColumnIJ_SeqAIJ, 2131 MatFDColoringCreate_SeqAIJ, 2132 0, 2133 0, 2134 MatPermute_SeqAIJ, 2135 0, 2136 0, 2137 MatDestroy_SeqAIJ, 2138 MatView_SeqAIJ, 2139 MatGetPetscMaps_Petsc, 2140 0, 2141 0, 2142 0, 2143 0, 2144 0, 2145 0, 2146 0, 2147 0, 2148 MatSetColoring_SeqAIJ, 2149 MatSetValuesAdic_SeqAIJ, 2150 MatSetValuesAdifor_SeqAIJ, 2151 MatFDColoringApply_SeqAIJ}; 2152 2153 EXTERN_C_BEGIN 2154 #undef __FUNCT__ 2155 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2156 2157 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2158 { 2159 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2160 int i,nz,n; 2161 2162 PetscFunctionBegin; 2163 2164 nz = aij->maxnz; 2165 n = mat->n; 2166 for (i=0; i<nz; i++) { 2167 aij->j[i] = indices[i]; 2168 } 2169 aij->nz = nz; 2170 for (i=0; i<n; i++) { 2171 aij->ilen[i] = aij->imax[i]; 2172 } 2173 2174 PetscFunctionReturn(0); 2175 } 2176 EXTERN_C_END 2177 2178 #undef __FUNCT__ 2179 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2180 /*@ 2181 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2182 in the matrix. 2183 2184 Input Parameters: 2185 + mat - the SeqAIJ matrix 2186 - indices - the column indices 2187 2188 Level: advanced 2189 2190 Notes: 2191 This can be called if you have precomputed the nonzero structure of the 2192 matrix and want to provide it to the matrix object to improve the performance 2193 of the MatSetValues() operation. 2194 2195 You MUST have set the correct numbers of nonzeros per row in the call to 2196 MatCreateSeqAIJ(). 2197 2198 MUST be called before any calls to MatSetValues(); 2199 2200 The indices should start with zero, not one. 2201 2202 @*/ 2203 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2204 { 2205 int ierr,(*f)(Mat,int *); 2206 2207 PetscFunctionBegin; 2208 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2209 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2210 if (f) { 2211 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2212 } else { 2213 SETERRQ(1,"Wrong type of matrix to set column indices"); 2214 } 2215 PetscFunctionReturn(0); 2216 } 2217 2218 /* ----------------------------------------------------------------------------------------*/ 2219 2220 EXTERN_C_BEGIN 2221 #undef __FUNCT__ 2222 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2223 int MatStoreValues_SeqAIJ(Mat mat) 2224 { 2225 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2226 size_t nz = aij->i[mat->m],ierr; 2227 2228 PetscFunctionBegin; 2229 if (aij->nonew != 1) { 2230 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2231 } 2232 2233 /* allocate space for values if not already there */ 2234 if (!aij->saved_values) { 2235 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2236 } 2237 2238 /* copy values over */ 2239 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2240 PetscFunctionReturn(0); 2241 } 2242 EXTERN_C_END 2243 2244 #undef __FUNCT__ 2245 #define __FUNCT__ "MatStoreValues" 2246 /*@ 2247 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2248 example, reuse of the linear part of a Jacobian, while recomputing the 2249 nonlinear portion. 2250 2251 Collect on Mat 2252 2253 Input Parameters: 2254 . mat - the matrix (currently on AIJ matrices support this option) 2255 2256 Level: advanced 2257 2258 Common Usage, with SNESSolve(): 2259 $ Create Jacobian matrix 2260 $ Set linear terms into matrix 2261 $ Apply boundary conditions to matrix, at this time matrix must have 2262 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2263 $ boundary conditions again will not change the nonzero structure 2264 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2265 $ ierr = MatStoreValues(mat); 2266 $ Call SNESSetJacobian() with matrix 2267 $ In your Jacobian routine 2268 $ ierr = MatRetrieveValues(mat); 2269 $ Set nonlinear terms in matrix 2270 2271 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2272 $ // build linear portion of Jacobian 2273 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2274 $ ierr = MatStoreValues(mat); 2275 $ loop over nonlinear iterations 2276 $ ierr = MatRetrieveValues(mat); 2277 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2278 $ // call MatAssemblyBegin/End() on matrix 2279 $ Solve linear system with Jacobian 2280 $ endloop 2281 2282 Notes: 2283 Matrix must already be assemblied before calling this routine 2284 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2285 calling this routine. 2286 2287 .seealso: MatRetrieveValues() 2288 2289 @*/ 2290 int MatStoreValues(Mat mat) 2291 { 2292 int ierr,(*f)(Mat); 2293 2294 PetscFunctionBegin; 2295 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2296 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2297 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2298 2299 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2300 if (f) { 2301 ierr = (*f)(mat);CHKERRQ(ierr); 2302 } else { 2303 SETERRQ(1,"Wrong type of matrix to store values"); 2304 } 2305 PetscFunctionReturn(0); 2306 } 2307 2308 EXTERN_C_BEGIN 2309 #undef __FUNCT__ 2310 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2311 int MatRetrieveValues_SeqAIJ(Mat mat) 2312 { 2313 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2314 int nz = aij->i[mat->m],ierr; 2315 2316 PetscFunctionBegin; 2317 if (aij->nonew != 1) { 2318 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2319 } 2320 if (!aij->saved_values) { 2321 SETERRQ(1,"Must call MatStoreValues(A);first"); 2322 } 2323 2324 /* copy values over */ 2325 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 EXTERN_C_END 2329 2330 #undef __FUNCT__ 2331 #define __FUNCT__ "MatRetrieveValues" 2332 /*@ 2333 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2334 example, reuse of the linear part of a Jacobian, while recomputing the 2335 nonlinear portion. 2336 2337 Collect on Mat 2338 2339 Input Parameters: 2340 . mat - the matrix (currently on AIJ matrices support this option) 2341 2342 Level: advanced 2343 2344 .seealso: MatStoreValues() 2345 2346 @*/ 2347 int MatRetrieveValues(Mat mat) 2348 { 2349 int ierr,(*f)(Mat); 2350 2351 PetscFunctionBegin; 2352 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2353 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2354 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2355 2356 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2357 if (f) { 2358 ierr = (*f)(mat);CHKERRQ(ierr); 2359 } else { 2360 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2361 } 2362 PetscFunctionReturn(0); 2363 } 2364 2365 /* 2366 This allows SeqAIJ matrices to be passed to the matlab engine 2367 */ 2368 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2369 #include "engine.h" /* Matlab include file */ 2370 #include "mex.h" /* Matlab include file */ 2371 EXTERN_C_BEGIN 2372 #undef __FUNCT__ 2373 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2374 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2375 { 2376 int ierr,i,*ai,*aj; 2377 Mat B = (Mat)obj; 2378 mxArray *mat; 2379 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2380 2381 PetscFunctionBegin; 2382 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2383 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2384 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2385 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2386 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2387 2388 /* Matlab indices start at 0 for sparse (what a surprise) */ 2389 2390 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2391 mxSetName(mat,obj->name); 2392 engPutArray((Engine *)mengine,mat); 2393 PetscFunctionReturn(0); 2394 } 2395 EXTERN_C_END 2396 2397 EXTERN_C_BEGIN 2398 #undef __FUNCT__ 2399 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2400 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2401 { 2402 int ierr,ii; 2403 Mat mat = (Mat)obj; 2404 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2405 mxArray *mmat; 2406 2407 PetscFunctionBegin; 2408 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2409 2410 mmat = engGetArray((Engine *)mengine,obj->name); 2411 2412 aij->nz = (mxGetJc(mmat))[mat->m]; 2413 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2414 aij->j = (int*)(aij->a + aij->nz); 2415 aij->i = aij->j + aij->nz; 2416 aij->singlemalloc = PETSC_TRUE; 2417 aij->freedata = PETSC_TRUE; 2418 2419 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2420 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2421 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2422 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2423 2424 for (ii=0; ii<mat->m; ii++) { 2425 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2426 } 2427 2428 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2429 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2430 2431 PetscFunctionReturn(0); 2432 } 2433 EXTERN_C_END 2434 #endif 2435 2436 /* --------------------------------------------------------------------------------*/ 2437 #undef __FUNCT__ 2438 #define __FUNCT__ "MatCreateSeqAIJ" 2439 /*@C 2440 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2441 (the default parallel PETSc format). For good matrix assembly performance 2442 the user should preallocate the matrix storage by setting the parameter nz 2443 (or the array nnz). By setting these parameters accurately, performance 2444 during matrix assembly can be increased by more than a factor of 50. 2445 2446 Collective on MPI_Comm 2447 2448 Input Parameters: 2449 + comm - MPI communicator, set to PETSC_COMM_SELF 2450 . m - number of rows 2451 . n - number of columns 2452 . nz - number of nonzeros per row (same for all rows) 2453 - nnz - array containing the number of nonzeros in the various rows 2454 (possibly different for each row) or PETSC_NULL 2455 2456 Output Parameter: 2457 . A - the matrix 2458 2459 Notes: 2460 The AIJ format (also called the Yale sparse matrix format or 2461 compressed row storage), is fully compatible with standard Fortran 77 2462 storage. That is, the stored row and column indices can begin at 2463 either one (as in Fortran) or zero. See the users' manual for details. 2464 2465 Specify the preallocated storage with either nz or nnz (not both). 2466 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2467 allocation. For large problems you MUST preallocate memory or you 2468 will get TERRIBLE performance, see the users' manual chapter on matrices. 2469 2470 By default, this format uses inodes (identical nodes) when possible, to 2471 improve numerical efficiency of matrix-vector products and solves. We 2472 search for consecutive rows with the same nonzero structure, thereby 2473 reusing matrix information to achieve increased efficiency. 2474 2475 Options Database Keys: 2476 + -mat_aij_no_inode - Do not use inodes 2477 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2478 - -mat_aij_oneindex - Internally use indexing starting at 1 2479 rather than 0. Note that when calling MatSetValues(), 2480 the user still MUST index entries starting at 0! 2481 2482 Level: intermediate 2483 2484 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2485 2486 @*/ 2487 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2488 { 2489 int ierr; 2490 2491 PetscFunctionBegin; 2492 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2493 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2494 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 #define SKIP_ALLOCATION -4 2499 2500 #undef __FUNCT__ 2501 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2502 /*@C 2503 MatSeqAIJSetPreallocation - For good matrix assembly performance 2504 the user should preallocate the matrix storage by setting the parameter nz 2505 (or the array nnz). By setting these parameters accurately, performance 2506 during matrix assembly can be increased by more than a factor of 50. 2507 2508 Collective on MPI_Comm 2509 2510 Input Parameters: 2511 + comm - MPI communicator, set to PETSC_COMM_SELF 2512 . m - number of rows 2513 . n - number of columns 2514 . nz - number of nonzeros per row (same for all rows) 2515 - nnz - array containing the number of nonzeros in the various rows 2516 (possibly different for each row) or PETSC_NULL 2517 2518 Output Parameter: 2519 . A - the matrix 2520 2521 Notes: 2522 The AIJ format (also called the Yale sparse matrix format or 2523 compressed row storage), is fully compatible with standard Fortran 77 2524 storage. That is, the stored row and column indices can begin at 2525 either one (as in Fortran) or zero. See the users' manual for details. 2526 2527 Specify the preallocated storage with either nz or nnz (not both). 2528 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2529 allocation. For large problems you MUST preallocate memory or you 2530 will get TERRIBLE performance, see the users' manual chapter on matrices. 2531 2532 By default, this format uses inodes (identical nodes) when possible, to 2533 improve numerical efficiency of matrix-vector products and solves. We 2534 search for consecutive rows with the same nonzero structure, thereby 2535 reusing matrix information to achieve increased efficiency. 2536 2537 Options Database Keys: 2538 + -mat_aij_no_inode - Do not use inodes 2539 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2540 - -mat_aij_oneindex - Internally use indexing starting at 1 2541 rather than 0. Note that when calling MatSetValues(), 2542 the user still MUST index entries starting at 0! 2543 2544 Level: intermediate 2545 2546 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2547 2548 @*/ 2549 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2550 { 2551 int ierr,(*f)(Mat,int,const int[]); 2552 2553 PetscFunctionBegin; 2554 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2555 if (f) { 2556 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2557 } 2558 PetscFunctionReturn(0); 2559 } 2560 2561 EXTERN_C_BEGIN 2562 #undef __FUNCT__ 2563 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2564 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2565 { 2566 Mat_SeqAIJ *b; 2567 size_t len = 0; 2568 PetscTruth skipallocation = PETSC_FALSE; 2569 int i,ierr; 2570 2571 PetscFunctionBegin; 2572 2573 if (nz == SKIP_ALLOCATION) { 2574 skipallocation = PETSC_TRUE; 2575 nz = 0; 2576 } 2577 2578 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2579 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2580 if (nnz) { 2581 for (i=0; i<B->m; i++) { 2582 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2583 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); 2584 } 2585 } 2586 2587 B->preallocated = PETSC_TRUE; 2588 b = (Mat_SeqAIJ*)B->data; 2589 2590 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2591 if (!nnz) { 2592 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2593 else if (nz <= 0) nz = 1; 2594 for (i=0; i<B->m; i++) b->imax[i] = nz; 2595 nz = nz*B->m; 2596 } else { 2597 nz = 0; 2598 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2599 } 2600 2601 if (!skipallocation) { 2602 /* allocate the matrix space */ 2603 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2604 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2605 b->j = (int*)(b->a + nz); 2606 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2607 b->i = b->j + nz; 2608 b->i[0] = 0; 2609 for (i=1; i<B->m+1; i++) { 2610 b->i[i] = b->i[i-1] + b->imax[i-1]; 2611 } 2612 b->singlemalloc = PETSC_TRUE; 2613 b->freedata = PETSC_TRUE; 2614 } else { 2615 b->freedata = PETSC_FALSE; 2616 } 2617 2618 /* b->ilen will count nonzeros in each row so far. */ 2619 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2620 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2621 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2622 2623 b->nz = 0; 2624 b->maxnz = nz; 2625 B->info.nz_unneeded = (double)b->maxnz; 2626 PetscFunctionReturn(0); 2627 } 2628 EXTERN_C_END 2629 2630 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2631 2632 EXTERN_C_BEGIN 2633 EXTERN int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2634 EXTERN int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2635 EXTERN int MatReorderForNonzeroDiagonal_SeqAIJ(Mat,PetscReal,IS,IS); 2636 EXTERN int MatAdjustForInodes_SeqAIJ(Mat,IS*,IS*); 2637 EXTERN int MatSeqAIJGetInodeSizes_SeqAIJ(Mat,int*,int*[],int*); 2638 EXTERN_C_END 2639 2640 EXTERN_C_BEGIN 2641 #undef __FUNCT__ 2642 #define __FUNCT__ "MatCreate_SeqAIJ" 2643 int MatCreate_SeqAIJ(Mat B) 2644 { 2645 Mat_SeqAIJ *b; 2646 int ierr,size; 2647 PetscTruth flg; 2648 2649 PetscFunctionBegin; 2650 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2651 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2652 2653 B->m = B->M = PetscMax(B->m,B->M); 2654 B->n = B->N = PetscMax(B->n,B->N); 2655 2656 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2657 B->data = (void*)b; 2658 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2659 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2660 B->factor = 0; 2661 B->lupivotthreshold = 1.0; 2662 B->mapping = 0; 2663 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2664 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2665 b->row = 0; 2666 b->col = 0; 2667 b->icol = 0; 2668 b->reallocs = 0; 2669 2670 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2671 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2672 2673 b->sorted = PETSC_FALSE; 2674 b->ignorezeroentries = PETSC_FALSE; 2675 b->roworiented = PETSC_TRUE; 2676 b->nonew = 0; 2677 b->diag = 0; 2678 b->solve_work = 0; 2679 B->spptr = 0; 2680 b->inode.use = PETSC_TRUE; 2681 b->inode.node_count = 0; 2682 b->inode.size = 0; 2683 b->inode.limit = 5; 2684 b->inode.max_limit = 5; 2685 b->saved_values = 0; 2686 b->idiag = 0; 2687 b->ssor = 0; 2688 b->keepzeroedrows = PETSC_FALSE; 2689 b->xtoy = 0; 2690 b->XtoY = 0; 2691 2692 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2693 2694 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2695 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2696 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2697 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2698 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2699 if (flg){ 2700 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2701 } 2702 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2703 "MatSeqAIJSetColumnIndices_SeqAIJ", 2704 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2706 "MatStoreValues_SeqAIJ", 2707 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2709 "MatRetrieveValues_SeqAIJ", 2710 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2711 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2712 "MatConvert_SeqAIJ_SeqSBAIJ", 2713 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2714 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2715 "MatConvert_SeqAIJ_SeqBAIJ", 2716 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2717 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsSymmetric_C", 2718 "MatIsSymmetric_SeqAIJ", 2719 MatIsSymmetric_SeqAIJ);CHKERRQ(ierr); 2720 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2721 "MatSeqAIJSetPreallocation_SeqAIJ", 2722 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2723 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2724 "MatReorderForNonzeroDiagonal_SeqAIJ", 2725 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2726 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2727 "MatAdjustForInodes_SeqAIJ", 2728 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2730 "MatSeqAIJGetInodeSizes_SeqAIJ", 2731 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2732 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2735 #endif 2736 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2737 PetscFunctionReturn(0); 2738 } 2739 EXTERN_C_END 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2743 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2744 { 2745 Mat C; 2746 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2747 int i,m = A->m,ierr; 2748 size_t len; 2749 2750 PetscFunctionBegin; 2751 *B = 0; 2752 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2753 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2754 c = (Mat_SeqAIJ*)C->data; 2755 2756 C->factor = A->factor; 2757 c->row = 0; 2758 c->col = 0; 2759 c->icol = 0; 2760 c->keepzeroedrows = a->keepzeroedrows; 2761 C->assembled = PETSC_TRUE; 2762 2763 C->M = A->m; 2764 C->N = A->n; 2765 2766 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2767 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2768 for (i=0; i<m; i++) { 2769 c->imax[i] = a->imax[i]; 2770 c->ilen[i] = a->ilen[i]; 2771 } 2772 2773 /* allocate the matrix space */ 2774 c->singlemalloc = PETSC_TRUE; 2775 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2776 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2777 c->j = (int*)(c->a + a->i[m] ); 2778 c->i = c->j + a->i[m]; 2779 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2780 if (m > 0) { 2781 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2782 if (cpvalues == MAT_COPY_VALUES) { 2783 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2784 } else { 2785 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2786 } 2787 } 2788 2789 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2790 c->sorted = a->sorted; 2791 c->roworiented = a->roworiented; 2792 c->nonew = a->nonew; 2793 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2794 c->saved_values = 0; 2795 c->idiag = 0; 2796 c->ssor = 0; 2797 c->ignorezeroentries = a->ignorezeroentries; 2798 c->freedata = PETSC_TRUE; 2799 2800 if (a->diag) { 2801 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2802 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2803 for (i=0; i<m; i++) { 2804 c->diag[i] = a->diag[i]; 2805 } 2806 } else c->diag = 0; 2807 c->inode.use = a->inode.use; 2808 c->inode.limit = a->inode.limit; 2809 c->inode.max_limit = a->inode.max_limit; 2810 if (a->inode.size){ 2811 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2812 c->inode.node_count = a->inode.node_count; 2813 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2814 } else { 2815 c->inode.size = 0; 2816 c->inode.node_count = 0; 2817 } 2818 c->nz = a->nz; 2819 c->maxnz = a->maxnz; 2820 c->solve_work = 0; 2821 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2822 C->preallocated = PETSC_TRUE; 2823 2824 *B = C; 2825 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2826 PetscFunctionReturn(0); 2827 } 2828 2829 EXTERN_C_BEGIN 2830 #undef __FUNCT__ 2831 #define __FUNCT__ "MatLoad_SeqAIJ" 2832 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2833 { 2834 Mat_SeqAIJ *a; 2835 Mat B; 2836 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2837 MPI_Comm comm; 2838 #if defined(PETSC_HAVE_MUMPS) 2839 PetscTruth flag; 2840 #endif 2841 2842 PetscFunctionBegin; 2843 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2844 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2845 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2846 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2847 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2848 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2849 M = header[1]; N = header[2]; nz = header[3]; 2850 2851 if (nz < 0) { 2852 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2853 } 2854 2855 /* read in row lengths */ 2856 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2857 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2858 2859 /* create our matrix */ 2860 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2861 ierr = MatSetType(B,type);CHKERRQ(ierr); 2862 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2863 a = (Mat_SeqAIJ*)B->data; 2864 2865 /* read in column indices and adjust for Fortran indexing*/ 2866 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2867 2868 /* read in nonzero values */ 2869 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2870 2871 /* set matrix "i" values */ 2872 a->i[0] = 0; 2873 for (i=1; i<= M; i++) { 2874 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2875 a->ilen[i-1] = rowlengths[i-1]; 2876 } 2877 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2878 2879 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2880 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2881 #if defined(PETSC_HAVE_MUMPS) 2882 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_mumps",&flag);CHKERRQ(ierr); 2883 if (flag) { ierr = MatUseMUMPS_MPIAIJ(B);CHKERRQ(ierr); } 2884 #endif 2885 *A = B; 2886 PetscFunctionReturn(0); 2887 } 2888 EXTERN_C_END 2889 2890 #undef __FUNCT__ 2891 #define __FUNCT__ "MatEqual_SeqAIJ" 2892 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2893 { 2894 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2895 int ierr; 2896 2897 PetscFunctionBegin; 2898 /* If the matrix dimensions are not equal,or no of nonzeros */ 2899 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2900 *flg = PETSC_FALSE; 2901 PetscFunctionReturn(0); 2902 } 2903 2904 /* if the a->i are the same */ 2905 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2906 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2907 2908 /* if a->j are the same */ 2909 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2910 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2911 2912 /* if a->a are the same */ 2913 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2914 2915 PetscFunctionReturn(0); 2916 2917 } 2918 2919 #undef __FUNCT__ 2920 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2921 /*@C 2922 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2923 provided by the user. 2924 2925 Coolective on MPI_Comm 2926 2927 Input Parameters: 2928 + comm - must be an MPI communicator of size 1 2929 . m - number of rows 2930 . n - number of columns 2931 . i - row indices 2932 . j - column indices 2933 - a - matrix values 2934 2935 Output Parameter: 2936 . mat - the matrix 2937 2938 Level: intermediate 2939 2940 Notes: 2941 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2942 once the matrix is destroyed 2943 2944 You cannot set new nonzero locations into this matrix, that will generate an error. 2945 2946 The i and j indices are 0 based 2947 2948 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2949 2950 @*/ 2951 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2952 { 2953 int ierr,ii; 2954 Mat_SeqAIJ *aij; 2955 2956 PetscFunctionBegin; 2957 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2958 aij = (Mat_SeqAIJ*)(*mat)->data; 2959 2960 if (i[0] != 0) { 2961 SETERRQ(1,"i (row indices) must start with 0"); 2962 } 2963 aij->i = i; 2964 aij->j = j; 2965 aij->a = a; 2966 aij->singlemalloc = PETSC_FALSE; 2967 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2968 aij->freedata = PETSC_FALSE; 2969 2970 for (ii=0; ii<m; ii++) { 2971 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2972 #if defined(PETSC_USE_BOPT_g) 2973 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]); 2974 #endif 2975 } 2976 #if defined(PETSC_USE_BOPT_g) 2977 for (ii=0; ii<aij->i[m]; ii++) { 2978 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2979 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2980 } 2981 #endif 2982 2983 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2984 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2985 PetscFunctionReturn(0); 2986 } 2987 2988 #undef __FUNCT__ 2989 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2990 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2991 { 2992 int ierr; 2993 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2994 2995 PetscFunctionBegin; 2996 if (coloring->ctype == IS_COLORING_LOCAL) { 2997 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2998 a->coloring = coloring; 2999 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3000 int i,*larray; 3001 ISColoring ocoloring; 3002 ISColoringValue *colors; 3003 3004 /* set coloring for diagonal portion */ 3005 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3006 for (i=0; i<A->n; i++) { 3007 larray[i] = i; 3008 } 3009 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3010 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3011 for (i=0; i<A->n; i++) { 3012 colors[i] = coloring->colors[larray[i]]; 3013 } 3014 ierr = PetscFree(larray);CHKERRQ(ierr); 3015 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3016 a->coloring = ocoloring; 3017 } 3018 PetscFunctionReturn(0); 3019 } 3020 3021 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3022 EXTERN_C_BEGIN 3023 #include "adic/ad_utils.h" 3024 EXTERN_C_END 3025 3026 #undef __FUNCT__ 3027 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3028 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3029 { 3030 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3031 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3032 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3033 ISColoringValue *color; 3034 3035 PetscFunctionBegin; 3036 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3037 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3038 color = a->coloring->colors; 3039 /* loop over rows */ 3040 for (i=0; i<m; i++) { 3041 nz = ii[i+1] - ii[i]; 3042 /* loop over columns putting computed value into matrix */ 3043 for (j=0; j<nz; j++) { 3044 *v++ = values[color[*jj++]]; 3045 } 3046 values += nlen; /* jump to next row of derivatives */ 3047 } 3048 PetscFunctionReturn(0); 3049 } 3050 3051 #else 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3055 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3056 { 3057 PetscFunctionBegin; 3058 SETERRQ(1,"PETSc installed without ADIC"); 3059 } 3060 3061 #endif 3062 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3065 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3066 { 3067 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3068 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3069 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3070 ISColoringValue *color; 3071 3072 PetscFunctionBegin; 3073 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3074 color = a->coloring->colors; 3075 /* loop over rows */ 3076 for (i=0; i<m; i++) { 3077 nz = ii[i+1] - ii[i]; 3078 /* loop over columns putting computed value into matrix */ 3079 for (j=0; j<nz; j++) { 3080 *v++ = values[color[*jj++]]; 3081 } 3082 values += nl; /* jump to next row of derivatives */ 3083 } 3084 PetscFunctionReturn(0); 3085 } 3086 3087