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