1 2 /*$Id: aij.c,v 1.382 2001/08/24 15:26:27 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 818 #undef __FUNCT__ 819 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 820 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 821 { 822 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 823 PetscScalar *x,*y; 824 int ierr,m = A->m,shift = a->indexshift; 825 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 826 PetscScalar *v,alpha; 827 int n,i,*idx; 828 #endif 829 830 PetscFunctionBegin; 831 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 832 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 833 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 834 y = y + shift; /* shift for Fortran start by 1 indexing */ 835 836 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 837 fortranmulttransposeaddaij_(&m,x,a->i,a->j+shift,a->a+shift,y); 838 #else 839 for (i=0; i<m; i++) { 840 idx = a->j + a->i[i] + shift; 841 v = a->a + a->i[i] + shift; 842 n = a->i[i+1] - a->i[i]; 843 alpha = x[i]; 844 while (n-->0) {y[*idx++] += alpha * *v++;} 845 } 846 #endif 847 PetscLogFlops(2*a->nz); 848 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 849 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 850 PetscFunctionReturn(0); 851 } 852 853 #undef __FUNCT__ 854 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 855 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 856 { 857 PetscScalar zero = 0.0;; 858 int ierr; 859 860 PetscFunctionBegin; 861 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 862 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 863 PetscFunctionReturn(0); 864 } 865 866 867 #undef __FUNCT__ 868 #define __FUNCT__ "MatMult_SeqAIJ" 869 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 870 { 871 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 872 PetscScalar *x,*y,*v,sum; 873 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 874 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 875 int n,i,jrow,j; 876 #endif 877 878 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 879 #pragma disjoint(*x,*y,*v) 880 #endif 881 882 PetscFunctionBegin; 883 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 884 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 885 x = x + shift; /* shift for Fortran start by 1 indexing */ 886 idx = a->j; 887 v = a->a; 888 ii = a->i; 889 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 890 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 891 #else 892 v += shift; /* shift for Fortran start by 1 indexing */ 893 idx += shift; 894 for (i=0; i<m; i++) { 895 jrow = ii[i]; 896 n = ii[i+1] - jrow; 897 sum = 0.0; 898 for (j=0; j<n; j++) { 899 sum += v[jrow]*x[idx[jrow]]; jrow++; 900 } 901 y[i] = sum; 902 } 903 #endif 904 PetscLogFlops(2*a->nz - m); 905 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 906 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 907 PetscFunctionReturn(0); 908 } 909 910 #undef __FUNCT__ 911 #define __FUNCT__ "MatMultAdd_SeqAIJ" 912 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 913 { 914 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 915 PetscScalar *x,*y,*z,*v,sum; 916 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 917 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 918 int n,i,jrow,j; 919 #endif 920 921 PetscFunctionBegin; 922 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 923 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 924 if (zz != yy) { 925 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 926 } else { 927 z = y; 928 } 929 x = x + shift; /* shift for Fortran start by 1 indexing */ 930 idx = a->j; 931 v = a->a; 932 ii = a->i; 933 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 934 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 935 #else 936 v += shift; /* shift for Fortran start by 1 indexing */ 937 idx += shift; 938 for (i=0; i<m; i++) { 939 jrow = ii[i]; 940 n = ii[i+1] - jrow; 941 sum = y[i]; 942 for (j=0; j<n; j++) { 943 sum += v[jrow]*x[idx[jrow]]; jrow++; 944 } 945 z[i] = sum; 946 } 947 #endif 948 PetscLogFlops(2*a->nz); 949 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 950 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 951 if (zz != yy) { 952 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 953 } 954 PetscFunctionReturn(0); 955 } 956 957 /* 958 Adds diagonal pointers to sparse matrix structure. 959 */ 960 #undef __FUNCT__ 961 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 962 int MatMarkDiagonal_SeqAIJ(Mat A) 963 { 964 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 965 int i,j,*diag,m = A->m,shift = a->indexshift,ierr; 966 967 PetscFunctionBegin; 968 if (a->diag) PetscFunctionReturn(0); 969 970 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 971 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 972 for (i=0; i<A->m; i++) { 973 diag[i] = a->i[i+1]; 974 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 975 if (a->j[j]+shift == i) { 976 diag[i] = j - shift; 977 break; 978 } 979 } 980 } 981 a->diag = diag; 982 PetscFunctionReturn(0); 983 } 984 985 /* 986 Checks for missing diagonals 987 */ 988 #undef __FUNCT__ 989 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 990 int MatMissingDiagonal_SeqAIJ(Mat A) 991 { 992 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 993 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 994 995 PetscFunctionBegin; 996 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 997 diag = a->diag; 998 for (i=0; i<A->m; i++) { 999 if (jj[diag[i]+shift] != i-shift) { 1000 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 1001 } 1002 } 1003 PetscFunctionReturn(0); 1004 } 1005 1006 #undef __FUNCT__ 1007 #define __FUNCT__ "MatRelax_SeqAIJ" 1008 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,Vec xx) 1009 { 1010 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1011 PetscScalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 1012 int ierr,*idx,*diag,n = A->n,m = A->m,i,shift = a->indexshift; 1013 1014 PetscFunctionBegin; 1015 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1016 if (xx != bb) { 1017 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1018 } else { 1019 b = x; 1020 } 1021 1022 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1023 diag = a->diag; 1024 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1025 if (flag == SOR_APPLY_UPPER) { 1026 /* apply (U + D/omega) to the vector */ 1027 bs = b + shift; 1028 for (i=0; i<m; i++) { 1029 d = fshift + a->a[diag[i] + shift]; 1030 n = a->i[i+1] - diag[i] - 1; 1031 PetscLogFlops(2*n-1); 1032 idx = a->j + diag[i] + (!shift); 1033 v = a->a + diag[i] + (!shift); 1034 sum = b[i]*d/omega; 1035 SPARSEDENSEDOT(sum,bs,v,idx,n); 1036 x[i] = sum; 1037 } 1038 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1039 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1040 PetscFunctionReturn(0); 1041 } 1042 1043 /* setup workspace for Eisenstat */ 1044 if (flag & SOR_EISENSTAT) { 1045 if (!a->idiag) { 1046 ierr = PetscMalloc(2*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1047 a->ssor = a->idiag + m; 1048 v = a->a; 1049 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1050 } 1051 t = a->ssor; 1052 idiag = a->idiag; 1053 } 1054 /* Let A = L + U + D; where L is lower trianglar, 1055 U is upper triangular, E is diagonal; This routine applies 1056 1057 (L + E)^{-1} A (U + E)^{-1} 1058 1059 to a vector efficiently using Eisenstat's trick. This is for 1060 the case of SSOR preconditioner, so E is D/omega where omega 1061 is the relaxation factor. 1062 */ 1063 1064 if (flag == SOR_APPLY_LOWER) { 1065 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1066 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1067 /* special case for omega = 1.0 saves flops and some integer ops */ 1068 PetscScalar *v2; 1069 1070 v2 = a->a; 1071 /* x = (E + U)^{-1} b */ 1072 for (i=m-1; i>=0; i--) { 1073 n = a->i[i+1] - diag[i] - 1; 1074 idx = a->j + diag[i] + 1; 1075 v = a->a + diag[i] + 1; 1076 sum = b[i]; 1077 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1078 x[i] = sum*idiag[i]; 1079 1080 /* t = b - (2*E - D)x */ 1081 t[i] = b[i] - (v2[diag[i]])*x[i]; 1082 } 1083 1084 /* t = (E + L)^{-1}t */ 1085 diag = a->diag; 1086 for (i=0; i<m; i++) { 1087 n = diag[i] - a->i[i]; 1088 idx = a->j + a->i[i]; 1089 v = a->a + a->i[i]; 1090 sum = t[i]; 1091 SPARSEDENSEMDOT(sum,t,v,idx,n); 1092 t[i] = sum*idiag[i]; 1093 1094 /* x = x + t */ 1095 x[i] += t[i]; 1096 } 1097 1098 PetscLogFlops(3*m-1 + 2*a->nz); 1099 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1100 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1101 PetscFunctionReturn(0); 1102 } else if (flag & SOR_EISENSTAT) { 1103 /* Let A = L + U + D; where L is lower trianglar, 1104 U is upper triangular, E is diagonal; This routine applies 1105 1106 (L + E)^{-1} A (U + E)^{-1} 1107 1108 to a vector efficiently using Eisenstat's trick. This is for 1109 the case of SSOR preconditioner, so E is D/omega where omega 1110 is the relaxation factor. 1111 */ 1112 scale = (2.0/omega) - 1.0; 1113 1114 /* x = (E + U)^{-1} b */ 1115 for (i=m-1; i>=0; i--) { 1116 d = fshift + a->a[diag[i] + shift]; 1117 n = a->i[i+1] - diag[i] - 1; 1118 idx = a->j + diag[i] + (!shift); 1119 v = a->a + diag[i] + (!shift); 1120 sum = b[i]; 1121 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1122 x[i] = omega*(sum/d); 1123 } 1124 1125 /* t = b - (2*E - D)x */ 1126 v = a->a; 1127 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1128 1129 /* t = (E + L)^{-1}t */ 1130 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1131 diag = a->diag; 1132 for (i=0; i<m; i++) { 1133 d = fshift + a->a[diag[i]+shift]; 1134 n = diag[i] - a->i[i]; 1135 idx = a->j + a->i[i] + shift; 1136 v = a->a + a->i[i] + shift; 1137 sum = t[i]; 1138 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1139 t[i] = omega*(sum/d); 1140 /* x = x + t */ 1141 x[i] += t[i]; 1142 } 1143 1144 PetscLogFlops(6*m-1 + 2*a->nz); 1145 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1146 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1147 PetscFunctionReturn(0); 1148 } 1149 if (flag & SOR_ZERO_INITIAL_GUESS) { 1150 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1151 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1152 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1153 #else 1154 for (i=0; i<m; i++) { 1155 d = fshift + a->a[diag[i]+shift]; 1156 n = diag[i] - a->i[i]; 1157 PetscLogFlops(2*n-1); 1158 idx = a->j + a->i[i] + shift; 1159 v = a->a + a->i[i] + shift; 1160 sum = b[i]; 1161 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1162 x[i] = omega*(sum/d); 1163 } 1164 #endif 1165 xb = x; 1166 } else xb = b; 1167 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1168 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1169 for (i=0; i<m; i++) { 1170 x[i] *= a->a[diag[i]+shift]; 1171 } 1172 PetscLogFlops(m); 1173 } 1174 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1175 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1176 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,xb); 1177 #else 1178 for (i=m-1; i>=0; i--) { 1179 d = fshift + a->a[diag[i] + shift]; 1180 n = a->i[i+1] - diag[i] - 1; 1181 PetscLogFlops(2*n-1); 1182 idx = a->j + diag[i] + (!shift); 1183 v = a->a + diag[i] + (!shift); 1184 sum = xb[i]; 1185 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1186 x[i] = omega*(sum/d); 1187 } 1188 #endif 1189 } 1190 its--; 1191 } 1192 while (its--) { 1193 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1194 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1195 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1196 #else 1197 for (i=0; i<m; i++) { 1198 d = fshift + a->a[diag[i]+shift]; 1199 n = a->i[i+1] - a->i[i]; 1200 PetscLogFlops(2*n-1); 1201 idx = a->j + a->i[i] + shift; 1202 v = a->a + a->i[i] + shift; 1203 sum = b[i]; 1204 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1205 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1206 } 1207 #endif 1208 } 1209 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1210 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1211 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1212 #else 1213 for (i=m-1; i>=0; i--) { 1214 d = fshift + a->a[diag[i] + shift]; 1215 n = a->i[i+1] - a->i[i]; 1216 PetscLogFlops(2*n-1); 1217 idx = a->j + a->i[i] + shift; 1218 v = a->a + a->i[i] + shift; 1219 sum = b[i]; 1220 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1221 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1222 } 1223 #endif 1224 } 1225 } 1226 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1227 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1233 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1234 { 1235 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1236 1237 PetscFunctionBegin; 1238 info->rows_global = (double)A->m; 1239 info->columns_global = (double)A->n; 1240 info->rows_local = (double)A->m; 1241 info->columns_local = (double)A->n; 1242 info->block_size = 1.0; 1243 info->nz_allocated = (double)a->maxnz; 1244 info->nz_used = (double)a->nz; 1245 info->nz_unneeded = (double)(a->maxnz - a->nz); 1246 info->assemblies = (double)A->num_ass; 1247 info->mallocs = (double)a->reallocs; 1248 info->memory = A->mem; 1249 if (A->factor) { 1250 info->fill_ratio_given = A->info.fill_ratio_given; 1251 info->fill_ratio_needed = A->info.fill_ratio_needed; 1252 info->factor_mallocs = A->info.factor_mallocs; 1253 } else { 1254 info->fill_ratio_given = 0; 1255 info->fill_ratio_needed = 0; 1256 info->factor_mallocs = 0; 1257 } 1258 PetscFunctionReturn(0); 1259 } 1260 1261 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1262 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1263 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 1264 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1265 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1266 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1267 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1268 1269 #undef __FUNCT__ 1270 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1271 int MatZeroRows_SeqAIJ(Mat A,IS is,PetscScalar *diag) 1272 { 1273 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1274 int i,ierr,N,*rows,m = A->m - 1,shift = a->indexshift; 1275 1276 PetscFunctionBegin; 1277 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1278 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1279 if (a->keepzeroedrows) { 1280 for (i=0; i<N; i++) { 1281 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1282 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1283 } 1284 if (diag) { 1285 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1286 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1287 for (i=0; i<N; i++) { 1288 a->a[a->diag[rows[i]]] = *diag; 1289 } 1290 } 1291 } else { 1292 if (diag) { 1293 for (i=0; i<N; i++) { 1294 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1295 if (a->ilen[rows[i]] > 0) { 1296 a->ilen[rows[i]] = 1; 1297 a->a[a->i[rows[i]]+shift] = *diag; 1298 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1299 } else { /* in case row was completely empty */ 1300 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1301 } 1302 } 1303 } else { 1304 for (i=0; i<N; i++) { 1305 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1306 a->ilen[rows[i]] = 0; 1307 } 1308 } 1309 } 1310 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1311 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1312 PetscFunctionReturn(0); 1313 } 1314 1315 #undef __FUNCT__ 1316 #define __FUNCT__ "MatGetOwnershipRange_SeqAIJ" 1317 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1318 { 1319 PetscFunctionBegin; 1320 if (m) *m = 0; 1321 if (n) *n = A->m; 1322 PetscFunctionReturn(0); 1323 } 1324 1325 #undef __FUNCT__ 1326 #define __FUNCT__ "MatGetRow_SeqAIJ" 1327 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1328 { 1329 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1330 int *itmp,i,shift = a->indexshift,ierr; 1331 1332 PetscFunctionBegin; 1333 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1334 1335 *nz = a->i[row+1] - a->i[row]; 1336 if (v) *v = a->a + a->i[row] + shift; 1337 if (idx) { 1338 itmp = a->j + a->i[row] + shift; 1339 if (*nz && shift) { 1340 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 1341 for (i=0; i<(*nz); i++) {(*idx)[i] = itmp[i] + shift;} 1342 } else if (*nz) { 1343 *idx = itmp; 1344 } 1345 else *idx = 0; 1346 } 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1352 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1353 { 1354 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1355 int ierr; 1356 1357 PetscFunctionBegin; 1358 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1359 PetscFunctionReturn(0); 1360 } 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "MatNorm_SeqAIJ" 1364 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *norm) 1365 { 1366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1367 PetscScalar *v = a->a; 1368 PetscReal sum = 0.0; 1369 int i,j,shift = a->indexshift,ierr; 1370 1371 PetscFunctionBegin; 1372 if (type == NORM_FROBENIUS) { 1373 for (i=0; i<a->nz; i++) { 1374 #if defined(PETSC_USE_COMPLEX) 1375 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1376 #else 1377 sum += (*v)*(*v); v++; 1378 #endif 1379 } 1380 *norm = sqrt(sum); 1381 } else if (type == NORM_1) { 1382 PetscReal *tmp; 1383 int *jj = a->j; 1384 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1385 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1386 *norm = 0.0; 1387 for (j=0; j<a->nz; j++) { 1388 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1389 } 1390 for (j=0; j<A->n; j++) { 1391 if (tmp[j] > *norm) *norm = tmp[j]; 1392 } 1393 ierr = PetscFree(tmp);CHKERRQ(ierr); 1394 } else if (type == NORM_INFINITY) { 1395 *norm = 0.0; 1396 for (j=0; j<A->m; j++) { 1397 v = a->a + a->i[j] + shift; 1398 sum = 0.0; 1399 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1400 sum += PetscAbsScalar(*v); v++; 1401 } 1402 if (sum > *norm) *norm = sum; 1403 } 1404 } else { 1405 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1406 } 1407 PetscFunctionReturn(0); 1408 } 1409 1410 #undef __FUNCT__ 1411 #define __FUNCT__ "MatTranspose_SeqAIJ" 1412 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1413 { 1414 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1415 Mat C; 1416 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1417 int shift = a->indexshift; 1418 PetscScalar *array = a->a; 1419 1420 PetscFunctionBegin; 1421 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1422 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1423 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1424 if (shift) { 1425 for (i=0; i<ai[m]-1; i++) aj[i] -= 1; 1426 } 1427 for (i=0; i<ai[m]+shift; i++) col[aj[i]] += 1; 1428 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1429 ierr = PetscFree(col);CHKERRQ(ierr); 1430 for (i=0; i<m; i++) { 1431 len = ai[i+1]-ai[i]; 1432 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1433 array += len; 1434 aj += len; 1435 } 1436 if (shift) { 1437 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1438 } 1439 1440 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1441 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1442 1443 if (B) { 1444 *B = C; 1445 } else { 1446 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1447 } 1448 PetscFunctionReturn(0); 1449 } 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1453 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1454 { 1455 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1456 PetscScalar *l,*r,x,*v; 1457 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj,shift = a->indexshift; 1458 1459 PetscFunctionBegin; 1460 if (ll) { 1461 /* The local size is used so that VecMPI can be passed to this routine 1462 by MatDiagonalScale_MPIAIJ */ 1463 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1464 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1465 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1466 v = a->a; 1467 for (i=0; i<m; i++) { 1468 x = l[i]; 1469 M = a->i[i+1] - a->i[i]; 1470 for (j=0; j<M; j++) { (*v++) *= x;} 1471 } 1472 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1473 PetscLogFlops(nz); 1474 } 1475 if (rr) { 1476 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1477 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1478 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1479 v = a->a; jj = a->j; 1480 for (i=0; i<nz; i++) { 1481 (*v++) *= r[*jj++ + shift]; 1482 } 1483 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1484 PetscLogFlops(nz); 1485 } 1486 PetscFunctionReturn(0); 1487 } 1488 1489 #undef __FUNCT__ 1490 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1491 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1492 { 1493 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1494 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1495 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1496 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1497 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1498 PetscScalar *a_new,*mat_a; 1499 Mat C; 1500 PetscTruth stride; 1501 1502 PetscFunctionBegin; 1503 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1504 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1505 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1506 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1507 1508 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1509 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1510 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1511 1512 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1513 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1514 if (stride && step == 1) { 1515 /* special case of contiguous rows */ 1516 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1517 starts = lens + nrows; 1518 /* loop over new rows determining lens and starting points */ 1519 for (i=0; i<nrows; i++) { 1520 kstart = ai[irow[i]]+shift; 1521 kend = kstart + ailen[irow[i]]; 1522 for (k=kstart; k<kend; k++) { 1523 if (aj[k]+shift >= first) { 1524 starts[i] = k; 1525 break; 1526 } 1527 } 1528 sum = 0; 1529 while (k < kend) { 1530 if (aj[k++]+shift >= first+ncols) break; 1531 sum++; 1532 } 1533 lens[i] = sum; 1534 } 1535 /* create submatrix */ 1536 if (scall == MAT_REUSE_MATRIX) { 1537 int n_cols,n_rows; 1538 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1539 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1540 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1541 C = *B; 1542 } else { 1543 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1544 } 1545 c = (Mat_SeqAIJ*)C->data; 1546 1547 /* loop over rows inserting into submatrix */ 1548 a_new = c->a; 1549 j_new = c->j; 1550 i_new = c->i; 1551 i_new[0] = -shift; 1552 for (i=0; i<nrows; i++) { 1553 ii = starts[i]; 1554 lensi = lens[i]; 1555 for (k=0; k<lensi; k++) { 1556 *j_new++ = aj[ii+k] - first; 1557 } 1558 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1559 a_new += lensi; 1560 i_new[i+1] = i_new[i] + lensi; 1561 c->ilen[i] = lensi; 1562 } 1563 ierr = PetscFree(lens);CHKERRQ(ierr); 1564 } else { 1565 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1566 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1567 ssmap = smap + shift; 1568 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1569 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1570 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1571 /* determine lens of each row */ 1572 for (i=0; i<nrows; i++) { 1573 kstart = ai[irow[i]]+shift; 1574 kend = kstart + a->ilen[irow[i]]; 1575 lens[i] = 0; 1576 for (k=kstart; k<kend; k++) { 1577 if (ssmap[aj[k]]) { 1578 lens[i]++; 1579 } 1580 } 1581 } 1582 /* Create and fill new matrix */ 1583 if (scall == MAT_REUSE_MATRIX) { 1584 PetscTruth equal; 1585 1586 c = (Mat_SeqAIJ *)((*B)->data); 1587 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1588 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1589 if (!equal) { 1590 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1591 } 1592 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1593 C = *B; 1594 } else { 1595 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1596 } 1597 c = (Mat_SeqAIJ *)(C->data); 1598 for (i=0; i<nrows; i++) { 1599 row = irow[i]; 1600 kstart = ai[row]+shift; 1601 kend = kstart + a->ilen[row]; 1602 mat_i = c->i[i]+shift; 1603 mat_j = c->j + mat_i; 1604 mat_a = c->a + mat_i; 1605 mat_ilen = c->ilen + i; 1606 for (k=kstart; k<kend; k++) { 1607 if ((tcol=ssmap[a->j[k]])) { 1608 *mat_j++ = tcol - (!shift); 1609 *mat_a++ = a->a[k]; 1610 (*mat_ilen)++; 1611 1612 } 1613 } 1614 } 1615 /* Free work space */ 1616 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1617 ierr = PetscFree(smap);CHKERRQ(ierr); 1618 ierr = PetscFree(lens);CHKERRQ(ierr); 1619 } 1620 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1621 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1622 1623 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1624 *B = C; 1625 PetscFunctionReturn(0); 1626 } 1627 1628 /* 1629 */ 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1632 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1633 { 1634 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1635 int ierr; 1636 Mat outA; 1637 PetscTruth row_identity,col_identity; 1638 1639 PetscFunctionBegin; 1640 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1641 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1642 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1643 if (!row_identity || !col_identity) { 1644 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1645 } 1646 1647 outA = inA; 1648 inA->factor = FACTOR_LU; 1649 a->row = row; 1650 a->col = col; 1651 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1652 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1653 1654 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1655 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1656 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1657 PetscLogObjectParent(inA,a->icol); 1658 1659 if (!a->solve_work) { /* this matrix may have been factored before */ 1660 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1661 } 1662 1663 if (!a->diag) { 1664 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1665 } 1666 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1667 PetscFunctionReturn(0); 1668 } 1669 1670 #include "petscblaslapack.h" 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatScale_SeqAIJ" 1673 int MatScale_SeqAIJ(PetscScalar *alpha,Mat inA) 1674 { 1675 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1676 int one = 1; 1677 1678 PetscFunctionBegin; 1679 BLscal_(&a->nz,alpha,a->a,&one); 1680 PetscLogFlops(a->nz); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1686 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1687 { 1688 int ierr,i; 1689 1690 PetscFunctionBegin; 1691 if (scall == MAT_INITIAL_MATRIX) { 1692 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1693 } 1694 1695 for (i=0; i<n; i++) { 1696 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1697 } 1698 PetscFunctionReturn(0); 1699 } 1700 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1703 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1704 { 1705 PetscFunctionBegin; 1706 *bs = 1; 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1712 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1713 { 1714 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1715 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1716 int start,end,*ai,*aj; 1717 PetscBT table; 1718 1719 PetscFunctionBegin; 1720 shift = a->indexshift; 1721 m = A->m; 1722 ai = a->i; 1723 aj = a->j+shift; 1724 1725 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1726 1727 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1728 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1729 1730 for (i=0; i<is_max; i++) { 1731 /* Initialize the two local arrays */ 1732 isz = 0; 1733 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1734 1735 /* Extract the indices, assume there can be duplicate entries */ 1736 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1737 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1738 1739 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1740 for (j=0; j<n ; ++j){ 1741 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1742 } 1743 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1744 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1745 1746 k = 0; 1747 for (j=0; j<ov; j++){ /* for each overlap */ 1748 n = isz; 1749 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1750 row = nidx[k]; 1751 start = ai[row]; 1752 end = ai[row+1]; 1753 for (l = start; l<end ; l++){ 1754 val = aj[l] + shift; 1755 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1756 } 1757 } 1758 } 1759 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1760 } 1761 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1762 ierr = PetscFree(nidx);CHKERRQ(ierr); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 /* -------------------------------------------------------------- */ 1767 #undef __FUNCT__ 1768 #define __FUNCT__ "MatPermute_SeqAIJ" 1769 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1770 { 1771 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1772 PetscScalar *vwork; 1773 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1774 int *row,*col,*cnew,j,*lens; 1775 IS icolp,irowp; 1776 1777 PetscFunctionBegin; 1778 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1779 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1780 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1781 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1782 1783 /* determine lengths of permuted rows */ 1784 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1785 for (i=0; i<m; i++) { 1786 lens[row[i]] = a->i[i+1] - a->i[i]; 1787 } 1788 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1789 ierr = PetscFree(lens);CHKERRQ(ierr); 1790 1791 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1792 for (i=0; i<m; i++) { 1793 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1794 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1795 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1796 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1797 } 1798 ierr = PetscFree(cnew);CHKERRQ(ierr); 1799 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1800 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1801 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1802 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1803 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1804 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 #undef __FUNCT__ 1809 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1810 int MatPrintHelp_SeqAIJ(Mat A) 1811 { 1812 static PetscTruth called = PETSC_FALSE; 1813 MPI_Comm comm = A->comm; 1814 int ierr; 1815 1816 PetscFunctionBegin; 1817 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1818 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1819 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1820 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1821 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1822 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1823 #if defined(PETSC_HAVE_ESSL) 1824 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1825 #endif 1826 #if defined(PETSC_HAVE_LUSOL) 1827 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1828 #endif 1829 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1830 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1831 #endif 1832 PetscFunctionReturn(0); 1833 } 1834 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1835 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1836 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "MatCopy_SeqAIJ" 1839 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1840 { 1841 int ierr; 1842 PetscTruth flg; 1843 1844 PetscFunctionBegin; 1845 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1846 if (str == SAME_NONZERO_PATTERN && flg) { 1847 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1848 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1849 1850 if (a->i[A->m]+a->indexshift != b->i[B->m]+a->indexshift) { 1851 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1852 } 1853 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m]+a->indexshift)*sizeof(PetscScalar));CHKERRQ(ierr); 1854 } else { 1855 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1856 } 1857 PetscFunctionReturn(0); 1858 } 1859 1860 #undef __FUNCT__ 1861 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1862 int MatSetUpPreallocation_SeqAIJ(Mat A) 1863 { 1864 int ierr; 1865 1866 PetscFunctionBegin; 1867 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1868 PetscFunctionReturn(0); 1869 } 1870 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "MatGetArray_SeqAIJ" 1873 int MatGetArray_SeqAIJ(Mat A,PetscScalar **array) 1874 { 1875 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1876 PetscFunctionBegin; 1877 *array = a->a; 1878 PetscFunctionReturn(0); 1879 } 1880 1881 #undef __FUNCT__ 1882 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1883 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar **array) 1884 { 1885 PetscFunctionBegin; 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1891 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1892 { 1893 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1894 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1895 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1896 PetscScalar *vscale_array; 1897 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1898 Vec w1,w2,w3; 1899 void *fctx = coloring->fctx; 1900 PetscTruth flg; 1901 1902 PetscFunctionBegin; 1903 if (!coloring->w1) { 1904 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1905 PetscLogObjectParent(coloring,coloring->w1); 1906 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1907 PetscLogObjectParent(coloring,coloring->w2); 1908 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1909 PetscLogObjectParent(coloring,coloring->w3); 1910 } 1911 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1912 1913 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1914 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1915 if (flg) { 1916 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1917 } else { 1918 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1919 } 1920 1921 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1922 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1923 1924 /* 1925 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1926 coloring->F for the coarser grids from the finest 1927 */ 1928 if (coloring->F) { 1929 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1930 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1931 if (m1 != m2) { 1932 coloring->F = 0; 1933 } 1934 } 1935 1936 if (coloring->F) { 1937 w1 = coloring->F; 1938 coloring->F = 0; 1939 } else { 1940 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1941 } 1942 1943 /* 1944 Compute all the scale factors and share with other processors 1945 */ 1946 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1947 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1948 for (k=0; k<coloring->ncolors; k++) { 1949 /* 1950 Loop over each column associated with color adding the 1951 perturbation to the vector w3. 1952 */ 1953 for (l=0; l<coloring->ncolumns[k]; l++) { 1954 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1955 dx = xx[col]; 1956 if (dx == 0.0) dx = 1.0; 1957 #if !defined(PETSC_USE_COMPLEX) 1958 if (dx < umin && dx >= 0.0) dx = umin; 1959 else if (dx < 0.0 && dx > -umin) dx = -umin; 1960 #else 1961 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1962 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1963 #endif 1964 dx *= epsilon; 1965 vscale_array[col] = 1.0/dx; 1966 } 1967 } 1968 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1969 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1970 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1971 1972 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1973 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1974 1975 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1976 else vscaleforrow = coloring->columnsforrow; 1977 1978 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1979 /* 1980 Loop over each color 1981 */ 1982 for (k=0; k<coloring->ncolors; k++) { 1983 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1984 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1985 /* 1986 Loop over each column associated with color adding the 1987 perturbation to the vector w3. 1988 */ 1989 for (l=0; l<coloring->ncolumns[k]; l++) { 1990 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1991 dx = xx[col]; 1992 if (dx == 0.0) dx = 1.0; 1993 #if !defined(PETSC_USE_COMPLEX) 1994 if (dx < umin && dx >= 0.0) dx = umin; 1995 else if (dx < 0.0 && dx > -umin) dx = -umin; 1996 #else 1997 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1998 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1999 #endif 2000 dx *= epsilon; 2001 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2002 w3_array[col] += dx; 2003 } 2004 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2005 2006 /* 2007 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2008 */ 2009 2010 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2011 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2012 2013 /* 2014 Loop over rows of vector, putting results into Jacobian matrix 2015 */ 2016 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2017 for (l=0; l<coloring->nrows[k]; l++) { 2018 row = coloring->rows[k][l]; 2019 col = coloring->columnsforrow[k][l]; 2020 y[row] *= vscale_array[vscaleforrow[k][l]]; 2021 srow = row + start; 2022 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2023 } 2024 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2025 } 2026 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2027 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2028 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2029 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2030 PetscFunctionReturn(0); 2031 } 2032 2033 /* -------------------------------------------------------------------*/ 2034 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2035 MatGetRow_SeqAIJ, 2036 MatRestoreRow_SeqAIJ, 2037 MatMult_SeqAIJ, 2038 MatMultAdd_SeqAIJ, 2039 MatMultTranspose_SeqAIJ, 2040 MatMultTransposeAdd_SeqAIJ, 2041 MatSolve_SeqAIJ, 2042 MatSolveAdd_SeqAIJ, 2043 MatSolveTranspose_SeqAIJ, 2044 MatSolveTransposeAdd_SeqAIJ, 2045 MatLUFactor_SeqAIJ, 2046 0, 2047 MatRelax_SeqAIJ, 2048 MatTranspose_SeqAIJ, 2049 MatGetInfo_SeqAIJ, 2050 MatEqual_SeqAIJ, 2051 MatGetDiagonal_SeqAIJ, 2052 MatDiagonalScale_SeqAIJ, 2053 MatNorm_SeqAIJ, 2054 0, 2055 MatAssemblyEnd_SeqAIJ, 2056 MatCompress_SeqAIJ, 2057 MatSetOption_SeqAIJ, 2058 MatZeroEntries_SeqAIJ, 2059 MatZeroRows_SeqAIJ, 2060 MatLUFactorSymbolic_SeqAIJ, 2061 MatLUFactorNumeric_SeqAIJ, 2062 0, 2063 0, 2064 MatSetUpPreallocation_SeqAIJ, 2065 0, 2066 MatGetOwnershipRange_SeqAIJ, 2067 MatILUFactorSymbolic_SeqAIJ, 2068 0, 2069 MatGetArray_SeqAIJ, 2070 MatRestoreArray_SeqAIJ, 2071 MatDuplicate_SeqAIJ, 2072 0, 2073 0, 2074 MatILUFactor_SeqAIJ, 2075 0, 2076 0, 2077 MatGetSubMatrices_SeqAIJ, 2078 MatIncreaseOverlap_SeqAIJ, 2079 MatGetValues_SeqAIJ, 2080 MatCopy_SeqAIJ, 2081 MatPrintHelp_SeqAIJ, 2082 MatScale_SeqAIJ, 2083 0, 2084 0, 2085 MatILUDTFactor_SeqAIJ, 2086 MatGetBlockSize_SeqAIJ, 2087 MatGetRowIJ_SeqAIJ, 2088 MatRestoreRowIJ_SeqAIJ, 2089 MatGetColumnIJ_SeqAIJ, 2090 MatRestoreColumnIJ_SeqAIJ, 2091 MatFDColoringCreate_SeqAIJ, 2092 0, 2093 0, 2094 MatPermute_SeqAIJ, 2095 0, 2096 0, 2097 MatDestroy_SeqAIJ, 2098 MatView_SeqAIJ, 2099 MatGetPetscMaps_Petsc, 2100 0, 2101 0, 2102 0, 2103 0, 2104 0, 2105 0, 2106 0, 2107 0, 2108 MatSetColoring_SeqAIJ, 2109 MatSetValuesAdic_SeqAIJ, 2110 MatSetValuesAdifor_SeqAIJ, 2111 MatFDColoringApply_SeqAIJ}; 2112 2113 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 2114 EXTERN int MatUseEssl_SeqAIJ(Mat); 2115 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 2116 EXTERN int MatUseMatlab_SeqAIJ(Mat); 2117 EXTERN int MatUseDXML_SeqAIJ(Mat); 2118 2119 EXTERN_C_BEGIN 2120 #undef __FUNCT__ 2121 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2122 2123 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2124 { 2125 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2126 int i,nz,n; 2127 2128 PetscFunctionBegin; 2129 2130 nz = aij->maxnz; 2131 n = mat->n; 2132 for (i=0; i<nz; i++) { 2133 aij->j[i] = indices[i]; 2134 } 2135 aij->nz = nz; 2136 for (i=0; i<n; i++) { 2137 aij->ilen[i] = aij->imax[i]; 2138 } 2139 2140 PetscFunctionReturn(0); 2141 } 2142 EXTERN_C_END 2143 2144 #undef __FUNCT__ 2145 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2146 /*@ 2147 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2148 in the matrix. 2149 2150 Input Parameters: 2151 + mat - the SeqAIJ matrix 2152 - indices - the column indices 2153 2154 Level: advanced 2155 2156 Notes: 2157 This can be called if you have precomputed the nonzero structure of the 2158 matrix and want to provide it to the matrix object to improve the performance 2159 of the MatSetValues() operation. 2160 2161 You MUST have set the correct numbers of nonzeros per row in the call to 2162 MatCreateSeqAIJ(). 2163 2164 MUST be called before any calls to MatSetValues(); 2165 2166 The indices should start with zero, not one. 2167 2168 @*/ 2169 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2170 { 2171 int ierr,(*f)(Mat,int *); 2172 2173 PetscFunctionBegin; 2174 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2175 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)())&f);CHKERRQ(ierr); 2176 if (f) { 2177 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2178 } else { 2179 SETERRQ(1,"Wrong type of matrix to set column indices"); 2180 } 2181 PetscFunctionReturn(0); 2182 } 2183 2184 /* ----------------------------------------------------------------------------------------*/ 2185 2186 EXTERN_C_BEGIN 2187 #undef __FUNCT__ 2188 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2189 int MatStoreValues_SeqAIJ(Mat mat) 2190 { 2191 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2192 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2193 2194 PetscFunctionBegin; 2195 if (aij->nonew != 1) { 2196 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2197 } 2198 2199 /* allocate space for values if not already there */ 2200 if (!aij->saved_values) { 2201 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2202 } 2203 2204 /* copy values over */ 2205 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2206 PetscFunctionReturn(0); 2207 } 2208 EXTERN_C_END 2209 2210 #undef __FUNCT__ 2211 #define __FUNCT__ "MatStoreValues" 2212 /*@ 2213 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2214 example, reuse of the linear part of a Jacobian, while recomputing the 2215 nonlinear portion. 2216 2217 Collect on Mat 2218 2219 Input Parameters: 2220 . mat - the matrix (currently on AIJ matrices support this option) 2221 2222 Level: advanced 2223 2224 Common Usage, with SNESSolve(): 2225 $ Create Jacobian matrix 2226 $ Set linear terms into matrix 2227 $ Apply boundary conditions to matrix, at this time matrix must have 2228 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2229 $ boundary conditions again will not change the nonzero structure 2230 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2231 $ ierr = MatStoreValues(mat); 2232 $ Call SNESSetJacobian() with matrix 2233 $ In your Jacobian routine 2234 $ ierr = MatRetrieveValues(mat); 2235 $ Set nonlinear terms in matrix 2236 2237 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2238 $ // build linear portion of Jacobian 2239 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2240 $ ierr = MatStoreValues(mat); 2241 $ loop over nonlinear iterations 2242 $ ierr = MatRetrieveValues(mat); 2243 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2244 $ // call MatAssemblyBegin/End() on matrix 2245 $ Solve linear system with Jacobian 2246 $ endloop 2247 2248 Notes: 2249 Matrix must already be assemblied before calling this routine 2250 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2251 calling this routine. 2252 2253 .seealso: MatRetrieveValues() 2254 2255 @*/ 2256 int MatStoreValues(Mat mat) 2257 { 2258 int ierr,(*f)(Mat); 2259 2260 PetscFunctionBegin; 2261 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2262 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2263 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2264 2265 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)())&f);CHKERRQ(ierr); 2266 if (f) { 2267 ierr = (*f)(mat);CHKERRQ(ierr); 2268 } else { 2269 SETERRQ(1,"Wrong type of matrix to store values"); 2270 } 2271 PetscFunctionReturn(0); 2272 } 2273 2274 EXTERN_C_BEGIN 2275 #undef __FUNCT__ 2276 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2277 int MatRetrieveValues_SeqAIJ(Mat mat) 2278 { 2279 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2280 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2281 2282 PetscFunctionBegin; 2283 if (aij->nonew != 1) { 2284 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2285 } 2286 if (!aij->saved_values) { 2287 SETERRQ(1,"Must call MatStoreValues(A);first"); 2288 } 2289 2290 /* copy values over */ 2291 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2292 PetscFunctionReturn(0); 2293 } 2294 EXTERN_C_END 2295 2296 #undef __FUNCT__ 2297 #define __FUNCT__ "MatRetrieveValues" 2298 /*@ 2299 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2300 example, reuse of the linear part of a Jacobian, while recomputing the 2301 nonlinear portion. 2302 2303 Collect on Mat 2304 2305 Input Parameters: 2306 . mat - the matrix (currently on AIJ matrices support this option) 2307 2308 Level: advanced 2309 2310 .seealso: MatStoreValues() 2311 2312 @*/ 2313 int MatRetrieveValues(Mat mat) 2314 { 2315 int ierr,(*f)(Mat); 2316 2317 PetscFunctionBegin; 2318 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2319 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2320 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2321 2322 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)())&f);CHKERRQ(ierr); 2323 if (f) { 2324 ierr = (*f)(mat);CHKERRQ(ierr); 2325 } else { 2326 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2327 } 2328 PetscFunctionReturn(0); 2329 } 2330 2331 /* 2332 This allows SeqAIJ matrices to be passed to the matlab engine 2333 */ 2334 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2335 #include "engine.h" /* Matlab include file */ 2336 #include "mex.h" /* Matlab include file */ 2337 EXTERN_C_BEGIN 2338 #undef __FUNCT__ 2339 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2340 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2341 { 2342 int ierr,i,*ai,*aj; 2343 Mat B = (Mat)obj; 2344 PetscScalar*array; 2345 mxArray *mat; 2346 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2347 2348 PetscFunctionBegin; 2349 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2350 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2351 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2352 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2353 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2354 2355 /* Matlab indices start at 0 for sparse (what a surprise) */ 2356 if (aij->indexshift) { 2357 for (i=0; i<B->m+1; i++) { 2358 ai[i]--; 2359 } 2360 for (i=0; i<aij->nz; i++) { 2361 aj[i]--; 2362 } 2363 } 2364 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2365 mxSetName(mat,obj->name); 2366 engPutArray((Engine *)engine,mat); 2367 PetscFunctionReturn(0); 2368 } 2369 EXTERN_C_END 2370 2371 EXTERN_C_BEGIN 2372 #undef __FUNCT__ 2373 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2374 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2375 { 2376 int ierr,ii; 2377 Mat mat = (Mat)obj; 2378 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2379 mxArray *mmat; 2380 2381 PetscFunctionBegin; 2382 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2383 aij->indexshift = 0; 2384 2385 mmat = engGetArray((Engine *)engine,obj->name); 2386 2387 aij->nz = (mxGetJc(mmat))[mat->m]; 2388 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2389 aij->j = (int*)(aij->a + aij->nz); 2390 aij->i = aij->j + aij->nz; 2391 aij->singlemalloc = PETSC_TRUE; 2392 aij->freedata = PETSC_TRUE; 2393 2394 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2395 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2396 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2397 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2398 2399 for (ii=0; ii<mat->m; ii++) { 2400 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2401 } 2402 2403 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2404 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2405 2406 PetscFunctionReturn(0); 2407 } 2408 EXTERN_C_END 2409 #endif 2410 2411 /* --------------------------------------------------------------------------------*/ 2412 2413 #undef __FUNCT__ 2414 #define __FUNCT__ "MatCreateSeqAIJ" 2415 /*@C 2416 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2417 (the default parallel PETSc format). For good matrix assembly performance 2418 the user should preallocate the matrix storage by setting the parameter nz 2419 (or the array nnz). By setting these parameters accurately, performance 2420 during matrix assembly can be increased by more than a factor of 50. 2421 2422 Collective on MPI_Comm 2423 2424 Input Parameters: 2425 + comm - MPI communicator, set to PETSC_COMM_SELF 2426 . m - number of rows 2427 . n - number of columns 2428 . nz - number of nonzeros per row (same for all rows) 2429 - nnz - array containing the number of nonzeros in the various rows 2430 (possibly different for each row) or PETSC_NULL 2431 2432 Output Parameter: 2433 . A - the matrix 2434 2435 Notes: 2436 The AIJ format (also called the Yale sparse matrix format or 2437 compressed row storage), is fully compatible with standard Fortran 77 2438 storage. That is, the stored row and column indices can begin at 2439 either one (as in Fortran) or zero. See the users' manual for details. 2440 2441 Specify the preallocated storage with either nz or nnz (not both). 2442 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2443 allocation. For large problems you MUST preallocate memory or you 2444 will get TERRIBLE performance, see the users' manual chapter on matrices. 2445 2446 By default, this format uses inodes (identical nodes) when possible, to 2447 improve numerical efficiency of matrix-vector products and solves. We 2448 search for consecutive rows with the same nonzero structure, thereby 2449 reusing matrix information to achieve increased efficiency. 2450 2451 Options Database Keys: 2452 + -mat_aij_no_inode - Do not use inodes 2453 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2454 - -mat_aij_oneindex - Internally use indexing starting at 1 2455 rather than 0. Note that when calling MatSetValues(), 2456 the user still MUST index entries starting at 0! 2457 2458 Level: intermediate 2459 2460 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2461 2462 @*/ 2463 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2464 { 2465 int ierr; 2466 2467 PetscFunctionBegin; 2468 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2469 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2470 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2471 PetscFunctionReturn(0); 2472 } 2473 2474 #undef __FUNCT__ 2475 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2476 /*@C 2477 MatSeqAIJSetPreallocation - For good matrix assembly performance 2478 the user should preallocate the matrix storage by setting the parameter nz 2479 (or the array nnz). By setting these parameters accurately, performance 2480 during matrix assembly can be increased by more than a factor of 50. 2481 2482 Collective on MPI_Comm 2483 2484 Input Parameters: 2485 + comm - MPI communicator, set to PETSC_COMM_SELF 2486 . m - number of rows 2487 . n - number of columns 2488 . nz - number of nonzeros per row (same for all rows) 2489 - nnz - array containing the number of nonzeros in the various rows 2490 (possibly different for each row) or PETSC_NULL 2491 2492 Output Parameter: 2493 . A - the matrix 2494 2495 Notes: 2496 The AIJ format (also called the Yale sparse matrix format or 2497 compressed row storage), is fully compatible with standard Fortran 77 2498 storage. That is, the stored row and column indices can begin at 2499 either one (as in Fortran) or zero. See the users' manual for details. 2500 2501 Specify the preallocated storage with either nz or nnz (not both). 2502 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2503 allocation. For large problems you MUST preallocate memory or you 2504 will get TERRIBLE performance, see the users' manual chapter on matrices. 2505 2506 By default, this format uses inodes (identical nodes) when possible, to 2507 improve numerical efficiency of matrix-vector products and solves. We 2508 search for consecutive rows with the same nonzero structure, thereby 2509 reusing matrix information to achieve increased efficiency. 2510 2511 Options Database Keys: 2512 + -mat_aij_no_inode - Do not use inodes 2513 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2514 - -mat_aij_oneindex - Internally use indexing starting at 1 2515 rather than 0. Note that when calling MatSetValues(), 2516 the user still MUST index entries starting at 0! 2517 2518 Level: intermediate 2519 2520 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2521 2522 @*/ 2523 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2524 { 2525 Mat_SeqAIJ *b; 2526 int i,len,ierr; 2527 PetscTruth flg2; 2528 2529 PetscFunctionBegin; 2530 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2531 if (!flg2) PetscFunctionReturn(0); 2532 2533 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2534 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2535 if (nnz) { 2536 for (i=0; i<B->m; i++) { 2537 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2538 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); 2539 } 2540 } 2541 2542 B->preallocated = PETSC_TRUE; 2543 b = (Mat_SeqAIJ*)B->data; 2544 2545 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2546 if (!nnz) { 2547 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2548 else if (nz <= 0) nz = 1; 2549 for (i=0; i<B->m; i++) b->imax[i] = nz; 2550 nz = nz*B->m; 2551 } else { 2552 nz = 0; 2553 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2554 } 2555 2556 /* allocate the matrix space */ 2557 len = nz*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2558 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2559 b->j = (int*)(b->a + nz); 2560 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2561 b->i = b->j + nz; 2562 b->singlemalloc = PETSC_TRUE; 2563 b->freedata = PETSC_TRUE; 2564 2565 b->i[0] = -b->indexshift; 2566 for (i=1; i<B->m+1; i++) { 2567 b->i[i] = b->i[i-1] + b->imax[i-1]; 2568 } 2569 2570 /* b->ilen will count nonzeros in each row so far. */ 2571 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2572 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2573 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2574 2575 b->nz = 0; 2576 b->maxnz = nz; 2577 B->info.nz_unneeded = (double)b->maxnz; 2578 PetscFunctionReturn(0); 2579 } 2580 2581 EXTERN_C_BEGIN 2582 #undef __FUNCT__ 2583 #define __FUNCT__ "MatCreate_SeqAIJ" 2584 int MatCreate_SeqAIJ(Mat B) 2585 { 2586 Mat_SeqAIJ *b; 2587 int ierr,size; 2588 PetscTruth flg; 2589 2590 PetscFunctionBegin; 2591 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2592 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2593 2594 B->m = B->M = PetscMax(B->m,B->M); 2595 B->n = B->N = PetscMax(B->n,B->N); 2596 2597 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2598 B->data = (void*)b; 2599 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2600 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2601 B->factor = 0; 2602 B->lupivotthreshold = 1.0; 2603 B->mapping = 0; 2604 ierr = PetscOptionsGetReal(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2605 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2606 b->row = 0; 2607 b->col = 0; 2608 b->icol = 0; 2609 b->indexshift = 0; 2610 b->reallocs = 0; 2611 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2612 if (flg) b->indexshift = -1; 2613 2614 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2615 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2616 2617 b->sorted = PETSC_FALSE; 2618 b->ignorezeroentries = PETSC_FALSE; 2619 b->roworiented = PETSC_TRUE; 2620 b->nonew = 0; 2621 b->diag = 0; 2622 b->solve_work = 0; 2623 b->spptr = 0; 2624 b->inode.use = PETSC_TRUE; 2625 b->inode.node_count = 0; 2626 b->inode.size = 0; 2627 b->inode.limit = 5; 2628 b->inode.max_limit = 5; 2629 b->saved_values = 0; 2630 b->idiag = 0; 2631 b->ssor = 0; 2632 b->keepzeroedrows = PETSC_FALSE; 2633 2634 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2635 2636 #if defined(PETSC_HAVE_SUPERLU) 2637 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2638 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2639 #endif 2640 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2641 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2642 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2643 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2644 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2645 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2646 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2647 if (flg) { 2648 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2649 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2650 } 2651 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2652 "MatSeqAIJSetColumnIndices_SeqAIJ", 2653 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2654 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2655 "MatStoreValues_SeqAIJ", 2656 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2657 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2658 "MatRetrieveValues_SeqAIJ", 2659 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2660 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2661 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2662 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2663 #endif 2664 PetscFunctionReturn(0); 2665 } 2666 EXTERN_C_END 2667 2668 #undef __FUNCT__ 2669 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2670 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2671 { 2672 Mat C; 2673 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2674 int i,len,m = A->m,shift = a->indexshift,ierr; 2675 2676 PetscFunctionBegin; 2677 *B = 0; 2678 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2679 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2680 c = (Mat_SeqAIJ*)C->data; 2681 2682 C->factor = A->factor; 2683 c->row = 0; 2684 c->col = 0; 2685 c->icol = 0; 2686 c->indexshift = shift; 2687 c->keepzeroedrows = a->keepzeroedrows; 2688 C->assembled = PETSC_TRUE; 2689 2690 C->M = A->m; 2691 C->N = A->n; 2692 2693 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2694 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2695 for (i=0; i<m; i++) { 2696 c->imax[i] = a->imax[i]; 2697 c->ilen[i] = a->ilen[i]; 2698 } 2699 2700 /* allocate the matrix space */ 2701 c->singlemalloc = PETSC_TRUE; 2702 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2703 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2704 c->j = (int*)(c->a + a->i[m] + shift); 2705 c->i = c->j + a->i[m] + shift; 2706 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2707 if (m > 0) { 2708 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2709 if (cpvalues == MAT_COPY_VALUES) { 2710 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2711 } else { 2712 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2713 } 2714 } 2715 2716 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2717 c->sorted = a->sorted; 2718 c->roworiented = a->roworiented; 2719 c->nonew = a->nonew; 2720 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2721 c->saved_values = 0; 2722 c->idiag = 0; 2723 c->ssor = 0; 2724 c->ignorezeroentries = a->ignorezeroentries; 2725 c->freedata = PETSC_TRUE; 2726 2727 if (a->diag) { 2728 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2729 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2730 for (i=0; i<m; i++) { 2731 c->diag[i] = a->diag[i]; 2732 } 2733 } else c->diag = 0; 2734 c->inode.use = a->inode.use; 2735 c->inode.limit = a->inode.limit; 2736 c->inode.max_limit = a->inode.max_limit; 2737 if (a->inode.size){ 2738 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2739 c->inode.node_count = a->inode.node_count; 2740 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2741 } else { 2742 c->inode.size = 0; 2743 c->inode.node_count = 0; 2744 } 2745 c->nz = a->nz; 2746 c->maxnz = a->maxnz; 2747 c->solve_work = 0; 2748 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2749 C->preallocated = PETSC_TRUE; 2750 2751 *B = C; 2752 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2753 PetscFunctionReturn(0); 2754 } 2755 2756 EXTERN_C_BEGIN 2757 #undef __FUNCT__ 2758 #define __FUNCT__ "MatLoad_SeqAIJ" 2759 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2760 { 2761 Mat_SeqAIJ *a; 2762 Mat B; 2763 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2764 MPI_Comm comm; 2765 2766 PetscFunctionBegin; 2767 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2768 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2769 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2770 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2771 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2772 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2773 M = header[1]; N = header[2]; nz = header[3]; 2774 2775 if (nz < 0) { 2776 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2777 } 2778 2779 /* read in row lengths */ 2780 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2781 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2782 2783 /* create our matrix */ 2784 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2785 B = *A; 2786 a = (Mat_SeqAIJ*)B->data; 2787 shift = a->indexshift; 2788 2789 /* read in column indices and adjust for Fortran indexing*/ 2790 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2791 if (shift) { 2792 for (i=0; i<nz; i++) { 2793 a->j[i] += 1; 2794 } 2795 } 2796 2797 /* read in nonzero values */ 2798 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2799 2800 /* set matrix "i" values */ 2801 a->i[0] = -shift; 2802 for (i=1; i<= M; i++) { 2803 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2804 a->ilen[i-1] = rowlengths[i-1]; 2805 } 2806 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2807 2808 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2809 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2810 PetscFunctionReturn(0); 2811 } 2812 EXTERN_C_END 2813 2814 #undef __FUNCT__ 2815 #define __FUNCT__ "MatEqual_SeqAIJ" 2816 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2817 { 2818 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2819 int ierr; 2820 PetscTruth flag; 2821 2822 PetscFunctionBegin; 2823 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2824 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2825 2826 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2827 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2828 *flg = PETSC_FALSE; 2829 PetscFunctionReturn(0); 2830 } 2831 2832 /* if the a->i are the same */ 2833 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2834 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2835 2836 /* if a->j are the same */ 2837 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2838 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2839 2840 /* if a->a are the same */ 2841 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2842 2843 PetscFunctionReturn(0); 2844 2845 } 2846 2847 #undef __FUNCT__ 2848 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2849 /*@C 2850 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2851 provided by the user. 2852 2853 Coolective on MPI_Comm 2854 2855 Input Parameters: 2856 + comm - must be an MPI communicator of size 1 2857 . m - number of rows 2858 . n - number of columns 2859 . i - row indices 2860 . j - column indices 2861 - a - matrix values 2862 2863 Output Parameter: 2864 . mat - the matrix 2865 2866 Level: intermediate 2867 2868 Notes: 2869 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2870 once the matrix is destroyed 2871 2872 You cannot set new nonzero locations into this matrix, that will generate an error. 2873 2874 The i and j indices can be either 0- or 1 based 2875 2876 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2877 2878 @*/ 2879 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2880 { 2881 int ierr,ii; 2882 Mat_SeqAIJ *aij; 2883 2884 PetscFunctionBegin; 2885 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2886 aij = (Mat_SeqAIJ*)(*mat)->data; 2887 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2888 2889 if (i[0] == 1) { 2890 aij->indexshift = -1; 2891 } else if (i[0]) { 2892 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2893 } 2894 aij->i = i; 2895 aij->j = j; 2896 aij->a = a; 2897 aij->singlemalloc = PETSC_FALSE; 2898 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2899 aij->freedata = PETSC_FALSE; 2900 2901 for (ii=0; ii<m; ii++) { 2902 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2903 #if defined(PETSC_USE_BOPT_g) 2904 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]); 2905 #endif 2906 } 2907 #if defined(PETSC_USE_BOPT_g) 2908 for (ii=0; ii<aij->i[m]; ii++) { 2909 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2910 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2911 } 2912 #endif 2913 2914 /* changes indices to start at 0 */ 2915 if (i[0]) { 2916 aij->indexshift = 0; 2917 for (ii=0; ii<m; ii++) { 2918 i[ii]--; 2919 } 2920 for (ii=0; ii<i[m]; ii++) { 2921 j[ii]--; 2922 } 2923 } 2924 2925 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2926 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2927 PetscFunctionReturn(0); 2928 } 2929 2930 #undef __FUNCT__ 2931 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2932 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2933 { 2934 int ierr; 2935 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2936 2937 PetscFunctionBegin; 2938 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2939 a->coloring = coloring; 2940 PetscFunctionReturn(0); 2941 } 2942 2943 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2944 EXTERN_C_BEGIN 2945 #include "adic_utils.h" 2946 EXTERN_C_END 2947 2948 #undef __FUNCT__ 2949 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2950 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2951 { 2952 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2953 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 2954 PetscScalar*v = a->a,*values; 2955 char *cadvalues = (char *)advalues; 2956 2957 PetscFunctionBegin; 2958 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2959 nlen = my_AD_GetDerivTypeSize(); 2960 color = a->coloring->colors; 2961 /* loop over rows */ 2962 for (i=0; i<m; i++) { 2963 nz = ii[i+1] - ii[i]; 2964 /* loop over columns putting computed value into matrix */ 2965 values = my_AD_GetGradArray(cadvalues); 2966 for (j=0; j<nz; j++) { 2967 *v++ = values[color[*jj++]]; 2968 } 2969 cadvalues += nlen; /* jump to next row of derivatives */ 2970 } 2971 PetscFunctionReturn(0); 2972 } 2973 2974 #else 2975 2976 #undef __FUNCT__ 2977 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2978 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2979 { 2980 PetscFunctionBegin; 2981 SETERRQ(1,"PETSc installed without ADIC"); 2982 } 2983 2984 #endif 2985 2986 #undef __FUNCT__ 2987 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 2988 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 2989 { 2990 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2991 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 2992 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 2993 2994 PetscFunctionBegin; 2995 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2996 color = a->coloring->colors; 2997 /* loop over rows */ 2998 for (i=0; i<m; i++) { 2999 nz = ii[i+1] - ii[i]; 3000 /* loop over columns putting computed value into matrix */ 3001 for (j=0; j<nz; j++) { 3002 *v++ = values[color[*jj++]]; 3003 } 3004 values += nl; /* jump to next row of derivatives */ 3005 } 3006 PetscFunctionReturn(0); 3007 } 3008 3009