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