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