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