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