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