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