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