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