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