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