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