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