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