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,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 671 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "MatCompress_SeqAIJ" 677 PetscErrorCode MatCompress_SeqAIJ(Mat A) 678 { 679 PetscFunctionBegin; 680 PetscFunctionReturn(0); 681 } 682 683 #undef __FUNCT__ 684 #define __FUNCT__ "MatSetOption_SeqAIJ" 685 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op) 686 { 687 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 688 PetscErrorCode ierr; 689 690 PetscFunctionBegin; 691 switch (op) { 692 case MAT_ROW_ORIENTED: 693 a->roworiented = PETSC_TRUE; 694 break; 695 case MAT_KEEP_ZEROED_ROWS: 696 a->keepzeroedrows = PETSC_TRUE; 697 break; 698 case MAT_COLUMN_ORIENTED: 699 a->roworiented = PETSC_FALSE; 700 break; 701 case MAT_COLUMNS_SORTED: 702 a->sorted = PETSC_TRUE; 703 break; 704 case MAT_COLUMNS_UNSORTED: 705 a->sorted = PETSC_FALSE; 706 break; 707 case MAT_NO_NEW_NONZERO_LOCATIONS: 708 a->nonew = 1; 709 break; 710 case MAT_NEW_NONZERO_LOCATION_ERR: 711 a->nonew = -1; 712 break; 713 case MAT_NEW_NONZERO_ALLOCATION_ERR: 714 a->nonew = -2; 715 break; 716 case MAT_YES_NEW_NONZERO_LOCATIONS: 717 a->nonew = 0; 718 break; 719 case MAT_IGNORE_ZERO_ENTRIES: 720 a->ignorezeroentries = PETSC_TRUE; 721 break; 722 case MAT_USE_COMPRESSEDROW: 723 a->compressedrow.use = PETSC_TRUE; 724 break; 725 case MAT_DO_NOT_USE_COMPRESSEDROW: 726 a->compressedrow.use = PETSC_FALSE; 727 break; 728 case MAT_ROWS_SORTED: 729 case MAT_ROWS_UNSORTED: 730 case MAT_YES_NEW_DIAGONALS: 731 case MAT_IGNORE_OFF_PROC_ENTRIES: 732 case MAT_USE_HASH_TABLE: 733 ierr = PetscLogInfo((A,"MatSetOption_SeqAIJ:Option ignored\n"));CHKERRQ(ierr); 734 break; 735 case MAT_NO_NEW_DIAGONALS: 736 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 737 default: 738 break; 739 } 740 ierr = MatSetOption_Inode(A,op);CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 746 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 747 { 748 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 749 PetscErrorCode ierr; 750 PetscInt i,j,n; 751 PetscScalar *x,zero = 0.0; 752 753 PetscFunctionBegin; 754 ierr = VecSet(v,zero);CHKERRQ(ierr); 755 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 756 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 757 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 758 for (i=0; i<A->m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 if (a->j[j] == i) { 761 x[i] = a->a[j]; 762 break; 763 } 764 } 765 } 766 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 767 PetscFunctionReturn(0); 768 } 769 770 #undef __FUNCT__ 771 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 772 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 773 { 774 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 775 PetscScalar *x,*y; 776 PetscErrorCode ierr; 777 PetscInt m = A->m; 778 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 779 PetscScalar *v,alpha; 780 PetscInt n,i,*idx,*ii,*ridx=PETSC_NULL; 781 Mat_CompressedRow cprow = a->compressedrow; 782 PetscTruth usecprow = cprow.use; 783 #endif 784 785 PetscFunctionBegin; 786 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 787 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 788 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 789 790 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 791 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 792 #else 793 if (usecprow){ 794 m = cprow.nrows; 795 ii = cprow.i; 796 ridx = cprow.rindex; 797 } else { 798 ii = a->i; 799 } 800 for (i=0; i<m; i++) { 801 idx = a->j + ii[i] ; 802 v = a->a + ii[i] ; 803 n = ii[i+1] - ii[i]; 804 if (usecprow){ 805 alpha = x[ridx[i]]; 806 } else { 807 alpha = x[i]; 808 } 809 while (n-->0) {y[*idx++] += alpha * *v++;} 810 } 811 #endif 812 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 813 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 814 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 815 PetscFunctionReturn(0); 816 } 817 818 #undef __FUNCT__ 819 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 820 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 821 { 822 PetscScalar zero = 0.0; 823 PetscErrorCode ierr; 824 825 PetscFunctionBegin; 826 ierr = VecSet(yy,zero);CHKERRQ(ierr); 827 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 832 #undef __FUNCT__ 833 #define __FUNCT__ "MatMult_SeqAIJ" 834 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 835 { 836 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 837 PetscScalar *x,*y,*aa; 838 PetscErrorCode ierr; 839 PetscInt m=A->m,*aj,*ii; 840 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 841 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 842 PetscScalar sum; 843 PetscTruth usecprow=a->compressedrow.use; 844 #endif 845 846 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 847 #pragma disjoint(*x,*y,*aa) 848 #endif 849 850 PetscFunctionBegin; 851 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 852 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 853 aj = a->j; 854 aa = a->a; 855 ii = a->i; 856 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 857 fortranmultaij_(&m,x,ii,aj,aa,y); 858 #else 859 if (usecprow){ /* use compressed row format */ 860 m = a->compressedrow.nrows; 861 ii = a->compressedrow.i; 862 ridx = a->compressedrow.rindex; 863 for (i=0; i<m; i++){ 864 n = ii[i+1] - ii[i]; 865 aj = a->j + ii[i]; 866 aa = a->a + ii[i]; 867 sum = 0.0; 868 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 869 y[*ridx++] = sum; 870 } 871 } else { /* do not use compressed row format */ 872 for (i=0; i<m; i++) { 873 jrow = ii[i]; 874 n = ii[i+1] - jrow; 875 sum = 0.0; 876 for (j=0; j<n; j++) { 877 sum += aa[jrow]*x[aj[jrow]]; jrow++; 878 } 879 y[i] = sum; 880 } 881 } 882 #endif 883 ierr = PetscLogFlops(2*a->nz - m);CHKERRQ(ierr); 884 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 885 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "MatMultAdd_SeqAIJ" 891 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 892 { 893 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 894 PetscScalar *x,*y,*z,*aa; 895 PetscErrorCode ierr; 896 PetscInt m = A->m,*aj,*ii; 897 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 898 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 899 PetscScalar sum; 900 PetscTruth usecprow=a->compressedrow.use; 901 #endif 902 903 PetscFunctionBegin; 904 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 905 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 906 if (zz != yy) { 907 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 908 } else { 909 z = y; 910 } 911 912 aj = a->j; 913 aa = a->a; 914 ii = a->i; 915 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 916 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 917 #else 918 if (usecprow){ /* use compressed row format */ 919 if (zz != yy){ 920 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 921 } 922 m = a->compressedrow.nrows; 923 ii = a->compressedrow.i; 924 ridx = a->compressedrow.rindex; 925 for (i=0; i<m; i++){ 926 n = ii[i+1] - ii[i]; 927 aj = a->j + ii[i]; 928 aa = a->a + ii[i]; 929 sum = y[*ridx]; 930 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 931 z[*ridx++] = sum; 932 } 933 } else { /* do not use compressed row format */ 934 for (i=0; i<m; i++) { 935 jrow = ii[i]; 936 n = ii[i+1] - jrow; 937 sum = y[i]; 938 for (j=0; j<n; j++) { 939 sum += aa[jrow]*x[aj[jrow]]; jrow++; 940 } 941 z[i] = sum; 942 } 943 } 944 #endif 945 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 946 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 947 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 948 if (zz != yy) { 949 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 950 } 951 PetscFunctionReturn(0); 952 } 953 954 /* 955 Adds diagonal pointers to sparse matrix structure. 956 */ 957 #undef __FUNCT__ 958 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 959 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 960 { 961 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 962 PetscErrorCode ierr; 963 PetscInt i,j,*diag,m = A->m; 964 965 PetscFunctionBegin; 966 if (a->diag) PetscFunctionReturn(0); 967 968 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&diag);CHKERRQ(ierr); 969 ierr = PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 970 for (i=0; i<A->m; i++) { 971 diag[i] = a->i[i+1]; 972 for (j=a->i[i]; j<a->i[i+1]; j++) { 973 if (a->j[j] == i) { 974 diag[i] = j; 975 break; 976 } 977 } 978 } 979 a->diag = diag; 980 PetscFunctionReturn(0); 981 } 982 983 /* 984 Checks for missing diagonals 985 */ 986 #undef __FUNCT__ 987 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 988 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A) 989 { 990 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 991 PetscErrorCode ierr; 992 PetscInt *diag,*jj = a->j,i; 993 994 PetscFunctionBegin; 995 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 996 diag = a->diag; 997 for (i=0; i<A->m; i++) { 998 if (jj[diag[i]] != i) { 999 SETERRQ1(PETSC_ERR_PLIB,"Matrix is missing diagonal number %D",i); 1000 } 1001 } 1002 PetscFunctionReturn(0); 1003 } 1004 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "MatRelax_SeqAIJ" 1007 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1008 { 1009 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1010 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1011 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1012 PetscErrorCode ierr; 1013 PetscInt n = A->n,m = A->m,i; 1014 const PetscInt *idx,*diag; 1015 1016 PetscFunctionBegin; 1017 its = its*lits; 1018 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1019 1020 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1021 diag = a->diag; 1022 if (!a->idiag) { 1023 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1024 a->ssor = a->idiag + m; 1025 mdiag = a->ssor + m; 1026 1027 v = a->a; 1028 1029 /* this is wrong when fshift omega changes each iteration */ 1030 if (omega == 1.0 && !fshift) { 1031 for (i=0; i<m; i++) { 1032 mdiag[i] = v[diag[i]]; 1033 a->idiag[i] = 1.0/v[diag[i]]; 1034 } 1035 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1036 } else { 1037 for (i=0; i<m; i++) { 1038 mdiag[i] = v[diag[i]]; 1039 a->idiag[i] = omega/(fshift + v[diag[i]]); 1040 } 1041 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1042 } 1043 } 1044 t = a->ssor; 1045 idiag = a->idiag; 1046 mdiag = a->idiag + 2*m; 1047 1048 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1049 if (xx != bb) { 1050 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1051 } else { 1052 b = x; 1053 } 1054 1055 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1056 xs = x; 1057 if (flag == SOR_APPLY_UPPER) { 1058 /* apply (U + D/omega) to the vector */ 1059 bs = b; 1060 for (i=0; i<m; i++) { 1061 d = fshift + a->a[diag[i]]; 1062 n = a->i[i+1] - diag[i] - 1; 1063 idx = a->j + diag[i] + 1; 1064 v = a->a + diag[i] + 1; 1065 sum = b[i]*d/omega; 1066 SPARSEDENSEDOT(sum,bs,v,idx,n); 1067 x[i] = sum; 1068 } 1069 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1070 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1071 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1072 PetscFunctionReturn(0); 1073 } 1074 1075 1076 /* Let A = L + U + D; where L is lower trianglar, 1077 U is upper triangular, E is diagonal; This routine applies 1078 1079 (L + E)^{-1} A (U + E)^{-1} 1080 1081 to a vector efficiently using Eisenstat's trick. This is for 1082 the case of SSOR preconditioner, so E is D/omega where omega 1083 is the relaxation factor. 1084 */ 1085 1086 if (flag == SOR_APPLY_LOWER) { 1087 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1088 } else if (flag & SOR_EISENSTAT) { 1089 /* Let A = L + U + D; where L is lower trianglar, 1090 U is upper triangular, E is diagonal; This routine applies 1091 1092 (L + E)^{-1} A (U + E)^{-1} 1093 1094 to a vector efficiently using Eisenstat's trick. This is for 1095 the case of SSOR preconditioner, so E is D/omega where omega 1096 is the relaxation factor. 1097 */ 1098 scale = (2.0/omega) - 1.0; 1099 1100 /* x = (E + U)^{-1} b */ 1101 for (i=m-1; i>=0; i--) { 1102 n = a->i[i+1] - diag[i] - 1; 1103 idx = a->j + diag[i] + 1; 1104 v = a->a + diag[i] + 1; 1105 sum = b[i]; 1106 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1107 x[i] = sum*idiag[i]; 1108 } 1109 1110 /* t = b - (2*E - D)x */ 1111 v = a->a; 1112 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1113 1114 /* t = (E + L)^{-1}t */ 1115 ts = t; 1116 diag = a->diag; 1117 for (i=0; i<m; i++) { 1118 n = diag[i] - a->i[i]; 1119 idx = a->j + a->i[i]; 1120 v = a->a + a->i[i]; 1121 sum = t[i]; 1122 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1123 t[i] = sum*idiag[i]; 1124 /* x = x + t */ 1125 x[i] += t[i]; 1126 } 1127 1128 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1129 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1130 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1131 PetscFunctionReturn(0); 1132 } 1133 if (flag & SOR_ZERO_INITIAL_GUESS) { 1134 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1135 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1136 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1137 #else 1138 for (i=0; i<m; i++) { 1139 n = diag[i] - a->i[i]; 1140 idx = a->j + a->i[i]; 1141 v = a->a + a->i[i]; 1142 sum = b[i]; 1143 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1144 x[i] = sum*idiag[i]; 1145 } 1146 #endif 1147 xb = x; 1148 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1149 } else xb = b; 1150 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1151 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1152 for (i=0; i<m; i++) { 1153 x[i] *= mdiag[i]; 1154 } 1155 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1156 } 1157 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1158 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1159 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1160 #else 1161 for (i=m-1; i>=0; i--) { 1162 n = a->i[i+1] - diag[i] - 1; 1163 idx = a->j + diag[i] + 1; 1164 v = a->a + diag[i] + 1; 1165 sum = xb[i]; 1166 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1167 x[i] = sum*idiag[i]; 1168 } 1169 #endif 1170 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1171 } 1172 its--; 1173 } 1174 while (its--) { 1175 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1176 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1177 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1178 #else 1179 for (i=0; i<m; i++) { 1180 d = fshift + a->a[diag[i]]; 1181 n = a->i[i+1] - a->i[i]; 1182 idx = a->j + a->i[i]; 1183 v = a->a + a->i[i]; 1184 sum = b[i]; 1185 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1186 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1187 } 1188 #endif 1189 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1190 } 1191 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1192 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1193 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1194 #else 1195 for (i=m-1; i>=0; i--) { 1196 d = fshift + a->a[diag[i]]; 1197 n = a->i[i+1] - a->i[i]; 1198 idx = a->j + a->i[i]; 1199 v = a->a + a->i[i]; 1200 sum = b[i]; 1201 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1202 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1203 } 1204 #endif 1205 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1206 } 1207 } 1208 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1209 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1215 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1216 { 1217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1218 1219 PetscFunctionBegin; 1220 info->rows_global = (double)A->m; 1221 info->columns_global = (double)A->n; 1222 info->rows_local = (double)A->m; 1223 info->columns_local = (double)A->n; 1224 info->block_size = 1.0; 1225 info->nz_allocated = (double)a->maxnz; 1226 info->nz_used = (double)a->nz; 1227 info->nz_unneeded = (double)(a->maxnz - a->nz); 1228 info->assemblies = (double)A->num_ass; 1229 info->mallocs = (double)a->reallocs; 1230 info->memory = A->mem; 1231 if (A->factor) { 1232 info->fill_ratio_given = A->info.fill_ratio_given; 1233 info->fill_ratio_needed = A->info.fill_ratio_needed; 1234 info->factor_mallocs = A->info.factor_mallocs; 1235 } else { 1236 info->fill_ratio_given = 0; 1237 info->fill_ratio_needed = 0; 1238 info->factor_mallocs = 0; 1239 } 1240 PetscFunctionReturn(0); 1241 } 1242 1243 #undef __FUNCT__ 1244 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1245 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1246 { 1247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1248 PetscInt i,m = A->m - 1; 1249 PetscErrorCode ierr; 1250 1251 PetscFunctionBegin; 1252 if (a->keepzeroedrows) { 1253 for (i=0; i<N; i++) { 1254 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1255 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1256 } 1257 if (diag != 0.0) { 1258 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1259 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1260 for (i=0; i<N; i++) { 1261 a->a[a->diag[rows[i]]] = diag; 1262 } 1263 } 1264 A->same_nonzero = PETSC_TRUE; 1265 } else { 1266 if (diag != 0.0) { 1267 for (i=0; i<N; i++) { 1268 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1269 if (a->ilen[rows[i]] > 0) { 1270 a->ilen[rows[i]] = 1; 1271 a->a[a->i[rows[i]]] = diag; 1272 a->j[a->i[rows[i]]] = rows[i]; 1273 } else { /* in case row was completely empty */ 1274 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1275 } 1276 } 1277 } else { 1278 for (i=0; i<N; i++) { 1279 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1280 a->ilen[rows[i]] = 0; 1281 } 1282 } 1283 A->same_nonzero = PETSC_FALSE; 1284 } 1285 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1286 PetscFunctionReturn(0); 1287 } 1288 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "MatGetRow_SeqAIJ" 1291 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1292 { 1293 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1294 PetscInt *itmp; 1295 1296 PetscFunctionBegin; 1297 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1298 1299 *nz = a->i[row+1] - a->i[row]; 1300 if (v) *v = a->a + a->i[row]; 1301 if (idx) { 1302 itmp = a->j + a->i[row]; 1303 if (*nz) { 1304 *idx = itmp; 1305 } 1306 else *idx = 0; 1307 } 1308 PetscFunctionReturn(0); 1309 } 1310 1311 /* remove this function? */ 1312 #undef __FUNCT__ 1313 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1314 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1315 { 1316 PetscFunctionBegin; 1317 PetscFunctionReturn(0); 1318 } 1319 1320 #undef __FUNCT__ 1321 #define __FUNCT__ "MatNorm_SeqAIJ" 1322 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1323 { 1324 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1325 PetscScalar *v = a->a; 1326 PetscReal sum = 0.0; 1327 PetscErrorCode ierr; 1328 PetscInt i,j; 1329 1330 PetscFunctionBegin; 1331 if (type == NORM_FROBENIUS) { 1332 for (i=0; i<a->nz; i++) { 1333 #if defined(PETSC_USE_COMPLEX) 1334 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1335 #else 1336 sum += (*v)*(*v); v++; 1337 #endif 1338 } 1339 *nrm = sqrt(sum); 1340 } else if (type == NORM_1) { 1341 PetscReal *tmp; 1342 PetscInt *jj = a->j; 1343 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1344 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1345 *nrm = 0.0; 1346 for (j=0; j<a->nz; j++) { 1347 tmp[*jj++] += PetscAbsScalar(*v); v++; 1348 } 1349 for (j=0; j<A->n; j++) { 1350 if (tmp[j] > *nrm) *nrm = tmp[j]; 1351 } 1352 ierr = PetscFree(tmp);CHKERRQ(ierr); 1353 } else if (type == NORM_INFINITY) { 1354 *nrm = 0.0; 1355 for (j=0; j<A->m; j++) { 1356 v = a->a + a->i[j]; 1357 sum = 0.0; 1358 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1359 sum += PetscAbsScalar(*v); v++; 1360 } 1361 if (sum > *nrm) *nrm = sum; 1362 } 1363 } else { 1364 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1365 } 1366 PetscFunctionReturn(0); 1367 } 1368 1369 #undef __FUNCT__ 1370 #define __FUNCT__ "MatTranspose_SeqAIJ" 1371 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1372 { 1373 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1374 Mat C; 1375 PetscErrorCode ierr; 1376 PetscInt i,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1377 PetscScalar *array = a->a; 1378 1379 PetscFunctionBegin; 1380 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1381 ierr = PetscMalloc((1+A->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1382 ierr = PetscMemzero(col,(1+A->n)*sizeof(PetscInt));CHKERRQ(ierr); 1383 1384 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1385 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1386 ierr = MatSetSizes(C,A->n,m,A->n,m);CHKERRQ(ierr); 1387 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1388 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1389 ierr = PetscFree(col);CHKERRQ(ierr); 1390 for (i=0; i<m; i++) { 1391 len = ai[i+1]-ai[i]; 1392 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1393 array += len; 1394 aj += len; 1395 } 1396 1397 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1398 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1399 1400 if (B) { 1401 *B = C; 1402 } else { 1403 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1404 } 1405 PetscFunctionReturn(0); 1406 } 1407 1408 EXTERN_C_BEGIN 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1411 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1412 { 1413 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1414 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1415 PetscErrorCode ierr; 1416 PetscInt ma,na,mb,nb, i; 1417 1418 PetscFunctionBegin; 1419 bij = (Mat_SeqAIJ *) B->data; 1420 1421 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1422 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1423 if (ma!=nb || na!=mb){ 1424 *f = PETSC_FALSE; 1425 PetscFunctionReturn(0); 1426 } 1427 aii = aij->i; bii = bij->i; 1428 adx = aij->j; bdx = bij->j; 1429 va = aij->a; vb = bij->a; 1430 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1431 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1432 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1433 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1434 1435 *f = PETSC_TRUE; 1436 for (i=0; i<ma; i++) { 1437 while (aptr[i]<aii[i+1]) { 1438 PetscInt idc,idr; 1439 PetscScalar vc,vr; 1440 /* column/row index/value */ 1441 idc = adx[aptr[i]]; 1442 idr = bdx[bptr[idc]]; 1443 vc = va[aptr[i]]; 1444 vr = vb[bptr[idc]]; 1445 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1446 *f = PETSC_FALSE; 1447 goto done; 1448 } else { 1449 aptr[i]++; 1450 if (B || i!=idc) bptr[idc]++; 1451 } 1452 } 1453 } 1454 done: 1455 ierr = PetscFree(aptr);CHKERRQ(ierr); 1456 if (B) { 1457 ierr = PetscFree(bptr);CHKERRQ(ierr); 1458 } 1459 PetscFunctionReturn(0); 1460 } 1461 EXTERN_C_END 1462 1463 #undef __FUNCT__ 1464 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1465 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1466 { 1467 PetscErrorCode ierr; 1468 PetscFunctionBegin; 1469 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1475 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1476 { 1477 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1478 PetscScalar *l,*r,x,*v; 1479 PetscErrorCode ierr; 1480 PetscInt i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1481 1482 PetscFunctionBegin; 1483 if (ll) { 1484 /* The local size is used so that VecMPI can be passed to this routine 1485 by MatDiagonalScale_MPIAIJ */ 1486 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1487 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1488 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1489 v = a->a; 1490 for (i=0; i<m; i++) { 1491 x = l[i]; 1492 M = a->i[i+1] - a->i[i]; 1493 for (j=0; j<M; j++) { (*v++) *= x;} 1494 } 1495 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1496 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1497 } 1498 if (rr) { 1499 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1500 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1501 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1502 v = a->a; jj = a->j; 1503 for (i=0; i<nz; i++) { 1504 (*v++) *= r[*jj++]; 1505 } 1506 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1507 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1508 } 1509 PetscFunctionReturn(0); 1510 } 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1514 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1515 { 1516 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1517 PetscErrorCode ierr; 1518 PetscInt *smap,i,k,kstart,kend,oldcols = A->n,*lens; 1519 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1520 PetscInt *irow,*icol,nrows,ncols; 1521 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1522 PetscScalar *a_new,*mat_a; 1523 Mat C; 1524 PetscTruth stride; 1525 1526 PetscFunctionBegin; 1527 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1528 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1529 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1530 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1531 1532 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1533 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1534 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1535 1536 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1537 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1538 if (stride && step == 1) { 1539 /* special case of contiguous rows */ 1540 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1541 starts = lens + nrows; 1542 /* loop over new rows determining lens and starting points */ 1543 for (i=0; i<nrows; i++) { 1544 kstart = ai[irow[i]]; 1545 kend = kstart + ailen[irow[i]]; 1546 for (k=kstart; k<kend; k++) { 1547 if (aj[k] >= first) { 1548 starts[i] = k; 1549 break; 1550 } 1551 } 1552 sum = 0; 1553 while (k < kend) { 1554 if (aj[k++] >= first+ncols) break; 1555 sum++; 1556 } 1557 lens[i] = sum; 1558 } 1559 /* create submatrix */ 1560 if (scall == MAT_REUSE_MATRIX) { 1561 PetscInt n_cols,n_rows; 1562 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1563 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1564 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1565 C = *B; 1566 } else { 1567 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1568 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1569 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1570 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1571 } 1572 c = (Mat_SeqAIJ*)C->data; 1573 1574 /* loop over rows inserting into submatrix */ 1575 a_new = c->a; 1576 j_new = c->j; 1577 i_new = c->i; 1578 1579 for (i=0; i<nrows; i++) { 1580 ii = starts[i]; 1581 lensi = lens[i]; 1582 for (k=0; k<lensi; k++) { 1583 *j_new++ = aj[ii+k] - first; 1584 } 1585 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1586 a_new += lensi; 1587 i_new[i+1] = i_new[i] + lensi; 1588 c->ilen[i] = lensi; 1589 } 1590 ierr = PetscFree(lens);CHKERRQ(ierr); 1591 } else { 1592 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1593 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1594 1595 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1596 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1597 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1598 /* determine lens of each row */ 1599 for (i=0; i<nrows; i++) { 1600 kstart = ai[irow[i]]; 1601 kend = kstart + a->ilen[irow[i]]; 1602 lens[i] = 0; 1603 for (k=kstart; k<kend; k++) { 1604 if (smap[aj[k]]) { 1605 lens[i]++; 1606 } 1607 } 1608 } 1609 /* Create and fill new matrix */ 1610 if (scall == MAT_REUSE_MATRIX) { 1611 PetscTruth equal; 1612 1613 c = (Mat_SeqAIJ *)((*B)->data); 1614 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1615 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1616 if (!equal) { 1617 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1618 } 1619 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(PetscInt));CHKERRQ(ierr); 1620 C = *B; 1621 } else { 1622 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1623 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1624 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1625 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1626 } 1627 c = (Mat_SeqAIJ *)(C->data); 1628 for (i=0; i<nrows; i++) { 1629 row = irow[i]; 1630 kstart = ai[row]; 1631 kend = kstart + a->ilen[row]; 1632 mat_i = c->i[i]; 1633 mat_j = c->j + mat_i; 1634 mat_a = c->a + mat_i; 1635 mat_ilen = c->ilen + i; 1636 for (k=kstart; k<kend; k++) { 1637 if ((tcol=smap[a->j[k]])) { 1638 *mat_j++ = tcol - 1; 1639 *mat_a++ = a->a[k]; 1640 (*mat_ilen)++; 1641 1642 } 1643 } 1644 } 1645 /* Free work space */ 1646 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1647 ierr = PetscFree(smap);CHKERRQ(ierr); 1648 ierr = PetscFree(lens);CHKERRQ(ierr); 1649 } 1650 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1651 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1652 1653 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1654 *B = C; 1655 PetscFunctionReturn(0); 1656 } 1657 1658 /* 1659 */ 1660 #undef __FUNCT__ 1661 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1662 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1663 { 1664 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1665 PetscErrorCode ierr; 1666 Mat outA; 1667 PetscTruth row_identity,col_identity; 1668 1669 PetscFunctionBegin; 1670 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1671 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1672 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1673 if (!row_identity || !col_identity) { 1674 SETERRQ(PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU"); 1675 } 1676 1677 outA = inA; 1678 inA->factor = FACTOR_LU; 1679 a->row = row; 1680 a->col = col; 1681 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1682 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1683 1684 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1685 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1686 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1687 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1688 1689 if (!a->solve_work) { /* this matrix may have been factored before */ 1690 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1691 } 1692 1693 if (!a->diag) { 1694 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1695 } 1696 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #include "petscblaslapack.h" 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "MatScale_SeqAIJ" 1703 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1704 { 1705 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1706 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1707 PetscScalar oalpha = alpha; 1708 PetscErrorCode ierr; 1709 1710 1711 PetscFunctionBegin; 1712 BLASscal_(&bnz,&oalpha,a->a,&one); 1713 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 1717 #undef __FUNCT__ 1718 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1719 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1720 { 1721 PetscErrorCode ierr; 1722 PetscInt i; 1723 1724 PetscFunctionBegin; 1725 if (scall == MAT_INITIAL_MATRIX) { 1726 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1727 } 1728 1729 for (i=0; i<n; i++) { 1730 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1731 } 1732 PetscFunctionReturn(0); 1733 } 1734 1735 #undef __FUNCT__ 1736 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1737 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1738 { 1739 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1740 PetscErrorCode ierr; 1741 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1742 PetscInt start,end,*ai,*aj; 1743 PetscBT table; 1744 1745 PetscFunctionBegin; 1746 m = A->m; 1747 ai = a->i; 1748 aj = a->j; 1749 1750 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1751 1752 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1753 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1754 1755 for (i=0; i<is_max; i++) { 1756 /* Initialize the two local arrays */ 1757 isz = 0; 1758 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1759 1760 /* Extract the indices, assume there can be duplicate entries */ 1761 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1762 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1763 1764 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1765 for (j=0; j<n ; ++j){ 1766 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1767 } 1768 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1769 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1770 1771 k = 0; 1772 for (j=0; j<ov; j++){ /* for each overlap */ 1773 n = isz; 1774 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1775 row = nidx[k]; 1776 start = ai[row]; 1777 end = ai[row+1]; 1778 for (l = start; l<end ; l++){ 1779 val = aj[l] ; 1780 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1781 } 1782 } 1783 } 1784 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1785 } 1786 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1787 ierr = PetscFree(nidx);CHKERRQ(ierr); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 /* -------------------------------------------------------------- */ 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "MatPermute_SeqAIJ" 1794 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1795 { 1796 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1797 PetscErrorCode ierr; 1798 PetscInt i,nz,m = A->m,n = A->n,*col; 1799 PetscInt *row,*cnew,j,*lens; 1800 IS icolp,irowp; 1801 PetscInt *cwork; 1802 PetscScalar *vwork; 1803 1804 PetscFunctionBegin; 1805 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1806 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1807 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1808 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1809 1810 /* determine lengths of permuted rows */ 1811 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1812 for (i=0; i<m; i++) { 1813 lens[row[i]] = a->i[i+1] - a->i[i]; 1814 } 1815 ierr = MatCreate(A->comm,B);CHKERRQ(ierr); 1816 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1817 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1818 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1819 ierr = PetscFree(lens);CHKERRQ(ierr); 1820 1821 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1822 for (i=0; i<m; i++) { 1823 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1824 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1825 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1826 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1827 } 1828 ierr = PetscFree(cnew);CHKERRQ(ierr); 1829 (*B)->assembled = PETSC_FALSE; 1830 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1831 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1832 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1833 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1834 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1835 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1836 PetscFunctionReturn(0); 1837 } 1838 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1841 PetscErrorCode MatPrintHelp_SeqAIJ(Mat A) 1842 { 1843 static PetscTruth called = PETSC_FALSE; 1844 MPI_Comm comm = A->comm; 1845 PetscErrorCode ierr; 1846 1847 PetscFunctionBegin; 1848 ierr = MatPrintHelp_Inode(A);CHKERRQ(ierr); 1849 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1850 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1851 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1852 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1853 ierr = (*PetscHelpPrintf)(comm," -mat_no_compressedrow: Do not use compressedrow\n");CHKERRQ(ierr); 1854 PetscFunctionReturn(0); 1855 } 1856 1857 #undef __FUNCT__ 1858 #define __FUNCT__ "MatCopy_SeqAIJ" 1859 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1860 { 1861 PetscErrorCode ierr; 1862 1863 PetscFunctionBegin; 1864 /* If the two matrices have the same copy implementation, use fast copy. */ 1865 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1866 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1867 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1868 1869 if (a->i[A->m] != b->i[B->m]) { 1870 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1871 } 1872 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1873 } else { 1874 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1875 } 1876 PetscFunctionReturn(0); 1877 } 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1881 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1882 { 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1887 PetscFunctionReturn(0); 1888 } 1889 1890 #undef __FUNCT__ 1891 #define __FUNCT__ "MatGetArray_SeqAIJ" 1892 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1893 { 1894 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1895 PetscFunctionBegin; 1896 *array = a->a; 1897 PetscFunctionReturn(0); 1898 } 1899 1900 #undef __FUNCT__ 1901 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1902 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1903 { 1904 PetscFunctionBegin; 1905 PetscFunctionReturn(0); 1906 } 1907 1908 #undef __FUNCT__ 1909 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1910 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1911 { 1912 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 1913 PetscErrorCode ierr; 1914 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1915 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1916 PetscScalar *vscale_array; 1917 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1918 Vec w1,w2,w3; 1919 void *fctx = coloring->fctx; 1920 PetscTruth flg; 1921 1922 PetscFunctionBegin; 1923 if (!coloring->w1) { 1924 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1925 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 1926 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1927 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 1928 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1929 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 1930 } 1931 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1932 1933 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1934 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1935 if (flg) { 1936 ierr = PetscLogInfo((coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"));CHKERRQ(ierr); 1937 } else { 1938 PetscTruth assembled; 1939 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 1940 if (assembled) { 1941 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1942 } 1943 } 1944 1945 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1946 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1947 1948 /* 1949 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1950 coloring->F for the coarser grids from the finest 1951 */ 1952 if (coloring->F) { 1953 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1954 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1955 if (m1 != m2) { 1956 coloring->F = 0; 1957 } 1958 } 1959 1960 if (coloring->F) { 1961 w1 = coloring->F; 1962 coloring->F = 0; 1963 } else { 1964 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1965 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1966 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1967 } 1968 1969 /* 1970 Compute all the scale factors and share with other processors 1971 */ 1972 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1973 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1974 for (k=0; k<coloring->ncolors; k++) { 1975 /* 1976 Loop over each column associated with color adding the 1977 perturbation to the vector w3. 1978 */ 1979 for (l=0; l<coloring->ncolumns[k]; l++) { 1980 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1981 dx = xx[col]; 1982 if (dx == 0.0) dx = 1.0; 1983 #if !defined(PETSC_USE_COMPLEX) 1984 if (dx < umin && dx >= 0.0) dx = umin; 1985 else if (dx < 0.0 && dx > -umin) dx = -umin; 1986 #else 1987 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1988 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1989 #endif 1990 dx *= epsilon; 1991 vscale_array[col] = 1.0/dx; 1992 } 1993 } 1994 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1995 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1996 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1997 1998 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1999 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2000 2001 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2002 else vscaleforrow = coloring->columnsforrow; 2003 2004 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2005 /* 2006 Loop over each color 2007 */ 2008 for (k=0; k<coloring->ncolors; k++) { 2009 coloring->currentcolor = k; 2010 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2011 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2012 /* 2013 Loop over each column associated with color adding the 2014 perturbation to the vector w3. 2015 */ 2016 for (l=0; l<coloring->ncolumns[k]; l++) { 2017 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2018 dx = xx[col]; 2019 if (dx == 0.0) dx = 1.0; 2020 #if !defined(PETSC_USE_COMPLEX) 2021 if (dx < umin && dx >= 0.0) dx = umin; 2022 else if (dx < 0.0 && dx > -umin) dx = -umin; 2023 #else 2024 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2025 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2026 #endif 2027 dx *= epsilon; 2028 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2029 w3_array[col] += dx; 2030 } 2031 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2032 2033 /* 2034 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2035 */ 2036 2037 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2038 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2039 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2040 ierr = VecAXPY(w2,mone,w1);CHKERRQ(ierr); 2041 2042 /* 2043 Loop over rows of vector, putting results into Jacobian matrix 2044 */ 2045 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2046 for (l=0; l<coloring->nrows[k]; l++) { 2047 row = coloring->rows[k][l]; 2048 col = coloring->columnsforrow[k][l]; 2049 y[row] *= vscale_array[vscaleforrow[k][l]]; 2050 srow = row + start; 2051 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2052 } 2053 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2054 } 2055 coloring->currentcolor = k; 2056 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2057 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2058 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2059 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063 #include "petscblaslapack.h" 2064 #undef __FUNCT__ 2065 #define __FUNCT__ "MatAXPY_SeqAIJ" 2066 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2067 { 2068 PetscErrorCode ierr; 2069 PetscInt i; 2070 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2071 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2072 2073 PetscFunctionBegin; 2074 if (str == SAME_NONZERO_PATTERN) { 2075 PetscScalar alpha = a; 2076 BLASaxpy_(&bnz,&alpha,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(Y,a,X,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(v,zero);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 #undef __FUNCT__ 2674 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2675 /*@C 2676 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2677 2678 Input Parameters: 2679 + B - the matrix 2680 . i - the indices into j for the start of each row (starts with zero) 2681 . j - the column indices for each row (starts with zero) these must be sorted for each row 2682 - v - optional values in the matrix 2683 2684 Contributed by: Lisandro Dalchin 2685 2686 Level: developer 2687 2688 .keywords: matrix, aij, compressed row, sparse, sequential 2689 2690 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2691 @*/ 2692 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2693 { 2694 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2695 PetscErrorCode ierr; 2696 2697 PetscFunctionBegin; 2698 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2699 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2700 if (f) { 2701 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2702 } 2703 PetscFunctionReturn(0); 2704 } 2705 2706 EXTERN_C_BEGIN 2707 #undef __FUNCT__ 2708 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2709 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt I[],const PetscInt J[],const PetscScalar v[]) 2710 { 2711 PetscInt i; 2712 PetscInt m,n; 2713 PetscInt nz; 2714 PetscInt *nnz, nz_max = 0; 2715 PetscScalar *values; 2716 PetscErrorCode ierr; 2717 2718 PetscFunctionBegin; 2719 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2720 2721 if (I[0]) { 2722 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "I[0] must be 0 it is %D", I[0]); 2723 } 2724 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2725 for(i = 0; i < m; i++) { 2726 nz = I[i+1]- I[i]; 2727 nz_max = PetscMax(nz_max, nz); 2728 if (nz < 0) { 2729 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2730 } 2731 nnz[i] = nz; 2732 } 2733 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2734 ierr = PetscFree(nnz);CHKERRQ(ierr); 2735 2736 if (v) { 2737 values = (PetscScalar*) v; 2738 } else { 2739 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2740 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2741 } 2742 2743 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 2744 2745 for(i = 0; i < m; i++) { 2746 nz = I[i+1] - I[i]; 2747 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+I[i], values + (v ? I[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2748 } 2749 2750 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2751 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2752 ierr = MatSetOption(B,MAT_COLUMNS_UNSORTED);CHKERRQ(ierr); 2753 2754 if (!v) { 2755 ierr = PetscFree(values);CHKERRQ(ierr); 2756 } 2757 PetscFunctionReturn(0); 2758 } 2759 EXTERN_C_END 2760 2761 /*MC 2762 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2763 based on compressed sparse row format. 2764 2765 Options Database Keys: 2766 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2767 2768 Level: beginner 2769 2770 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2771 M*/ 2772 2773 EXTERN_C_BEGIN 2774 #undef __FUNCT__ 2775 #define __FUNCT__ "MatCreate_SeqAIJ" 2776 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 2777 { 2778 Mat_SeqAIJ *b; 2779 PetscErrorCode ierr; 2780 PetscMPIInt size; 2781 2782 PetscFunctionBegin; 2783 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2784 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2785 2786 B->m = B->M = PetscMax(B->m,B->M); 2787 B->n = B->N = PetscMax(B->n,B->N); 2788 2789 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2790 B->data = (void*)b; 2791 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2792 B->factor = 0; 2793 B->mapping = 0; 2794 b->row = 0; 2795 b->col = 0; 2796 b->icol = 0; 2797 b->reallocs = 0; 2798 2799 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2800 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2801 2802 b->sorted = PETSC_FALSE; 2803 b->ignorezeroentries = PETSC_FALSE; 2804 b->roworiented = PETSC_TRUE; 2805 b->nonew = 0; 2806 b->diag = 0; 2807 b->solve_work = 0; 2808 B->spptr = 0; 2809 b->saved_values = 0; 2810 b->idiag = 0; 2811 b->ssor = 0; 2812 b->keepzeroedrows = PETSC_FALSE; 2813 b->xtoy = 0; 2814 b->XtoY = 0; 2815 b->compressedrow.use = PETSC_FALSE; 2816 b->compressedrow.nrows = B->m; 2817 b->compressedrow.i = PETSC_NULL; 2818 b->compressedrow.rindex = PETSC_NULL; 2819 b->compressedrow.checked = PETSC_FALSE; 2820 B->same_nonzero = PETSC_FALSE; 2821 2822 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2823 2824 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2825 "MatSeqAIJSetColumnIndices_SeqAIJ", 2826 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2827 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2828 "MatStoreValues_SeqAIJ", 2829 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2830 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2831 "MatRetrieveValues_SeqAIJ", 2832 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2833 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2834 "MatConvert_SeqAIJ_SeqSBAIJ", 2835 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2836 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2837 "MatConvert_SeqAIJ_SeqBAIJ", 2838 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2839 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2840 "MatIsTranspose_SeqAIJ", 2841 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2842 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2843 "MatSeqAIJSetPreallocation_SeqAIJ", 2844 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2845 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 2846 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 2847 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 2848 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2849 "MatReorderForNonzeroDiagonal_SeqAIJ", 2850 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2851 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 2852 PetscFunctionReturn(0); 2853 } 2854 EXTERN_C_END 2855 2856 #undef __FUNCT__ 2857 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2858 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2859 { 2860 Mat C; 2861 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2862 PetscErrorCode ierr; 2863 PetscInt i,m = A->m; 2864 2865 PetscFunctionBegin; 2866 *B = 0; 2867 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 2868 ierr = MatSetSizes(C,A->m,A->n,A->m,A->n);CHKERRQ(ierr); 2869 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2870 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2871 2872 c = (Mat_SeqAIJ*)C->data; 2873 2874 C->factor = A->factor; 2875 2876 c->row = 0; 2877 c->col = 0; 2878 c->icol = 0; 2879 c->reallocs = 0; 2880 2881 C->assembled = PETSC_TRUE; 2882 2883 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 2884 for (i=0; i<m; i++) { 2885 c->imax[i] = a->imax[i]; 2886 c->ilen[i] = a->ilen[i]; 2887 } 2888 2889 /* allocate the matrix space */ 2890 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 2891 c->singlemalloc = PETSC_TRUE; 2892 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2893 if (m > 0) { 2894 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2895 if (cpvalues == MAT_COPY_VALUES) { 2896 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2897 } else { 2898 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2899 } 2900 } 2901 2902 c->sorted = a->sorted; 2903 c->ignorezeroentries = a->ignorezeroentries; 2904 c->roworiented = a->roworiented; 2905 c->nonew = a->nonew; 2906 if (a->diag) { 2907 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2908 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2909 for (i=0; i<m; i++) { 2910 c->diag[i] = a->diag[i]; 2911 } 2912 } else c->diag = 0; 2913 c->solve_work = 0; 2914 c->saved_values = 0; 2915 c->idiag = 0; 2916 c->ssor = 0; 2917 c->keepzeroedrows = a->keepzeroedrows; 2918 c->freedata = PETSC_TRUE; 2919 c->xtoy = 0; 2920 c->XtoY = 0; 2921 2922 c->nz = a->nz; 2923 c->maxnz = a->maxnz; 2924 C->preallocated = PETSC_TRUE; 2925 2926 c->compressedrow.use = a->compressedrow.use; 2927 c->compressedrow.nrows = a->compressedrow.nrows; 2928 c->compressedrow.checked = a->compressedrow.checked; 2929 if ( a->compressedrow.checked && a->compressedrow.use){ 2930 i = a->compressedrow.nrows; 2931 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 2932 c->compressedrow.rindex = c->compressedrow.i + i + 1; 2933 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 2934 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 2935 } else { 2936 c->compressedrow.use = PETSC_FALSE; 2937 c->compressedrow.i = PETSC_NULL; 2938 c->compressedrow.rindex = PETSC_NULL; 2939 } 2940 C->same_nonzero = A->same_nonzero; 2941 2942 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 2943 2944 *B = C; 2945 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2946 PetscFunctionReturn(0); 2947 } 2948 2949 #undef __FUNCT__ 2950 #define __FUNCT__ "MatLoad_SeqAIJ" 2951 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 2952 { 2953 Mat_SeqAIJ *a; 2954 Mat B; 2955 PetscErrorCode ierr; 2956 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 2957 int fd; 2958 PetscMPIInt size; 2959 MPI_Comm comm; 2960 2961 PetscFunctionBegin; 2962 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2963 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2964 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2965 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2966 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2967 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2968 M = header[1]; N = header[2]; nz = header[3]; 2969 2970 if (nz < 0) { 2971 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2972 } 2973 2974 /* read in row lengths */ 2975 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2976 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2977 2978 /* check if sum of rowlengths is same as nz */ 2979 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 2980 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 2981 2982 /* create our matrix */ 2983 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2984 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2985 ierr = MatSetType(B,type);CHKERRQ(ierr); 2986 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 2987 a = (Mat_SeqAIJ*)B->data; 2988 2989 /* read in column indices and adjust for Fortran indexing*/ 2990 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2991 2992 /* read in nonzero values */ 2993 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2994 2995 /* set matrix "i" values */ 2996 a->i[0] = 0; 2997 for (i=1; i<= M; i++) { 2998 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2999 a->ilen[i-1] = rowlengths[i-1]; 3000 } 3001 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3002 3003 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3004 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3005 *A = B; 3006 PetscFunctionReturn(0); 3007 } 3008 3009 #undef __FUNCT__ 3010 #define __FUNCT__ "MatEqual_SeqAIJ" 3011 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3012 { 3013 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3014 PetscErrorCode ierr; 3015 3016 PetscFunctionBegin; 3017 /* If the matrix dimensions are not equal,or no of nonzeros */ 3018 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 3019 *flg = PETSC_FALSE; 3020 PetscFunctionReturn(0); 3021 } 3022 3023 /* if the a->i are the same */ 3024 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3025 if (!*flg) PetscFunctionReturn(0); 3026 3027 /* if a->j are the same */ 3028 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3029 if (!*flg) PetscFunctionReturn(0); 3030 3031 /* if a->a are the same */ 3032 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3033 3034 PetscFunctionReturn(0); 3035 3036 } 3037 3038 #undef __FUNCT__ 3039 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3040 /*@C 3041 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3042 provided by the user. 3043 3044 Coolective on MPI_Comm 3045 3046 Input Parameters: 3047 + comm - must be an MPI communicator of size 1 3048 . m - number of rows 3049 . n - number of columns 3050 . i - row indices 3051 . j - column indices 3052 - a - matrix values 3053 3054 Output Parameter: 3055 . mat - the matrix 3056 3057 Level: intermediate 3058 3059 Notes: 3060 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3061 once the matrix is destroyed 3062 3063 You cannot set new nonzero locations into this matrix, that will generate an error. 3064 3065 The i and j indices are 0 based 3066 3067 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 3068 3069 @*/ 3070 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3071 { 3072 PetscErrorCode ierr; 3073 PetscInt ii; 3074 Mat_SeqAIJ *aij; 3075 3076 PetscFunctionBegin; 3077 if (i[0]) { 3078 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3079 } 3080 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3081 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3082 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3083 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3084 aij = (Mat_SeqAIJ*)(*mat)->data; 3085 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3086 3087 aij->i = i; 3088 aij->j = j; 3089 aij->a = a; 3090 aij->singlemalloc = PETSC_FALSE; 3091 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3092 aij->freedata = PETSC_FALSE; 3093 3094 for (ii=0; ii<m; ii++) { 3095 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3096 #if defined(PETSC_USE_DEBUG) 3097 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]); 3098 #endif 3099 } 3100 #if defined(PETSC_USE_DEBUG) 3101 for (ii=0; ii<aij->i[m]; ii++) { 3102 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3103 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3104 } 3105 #endif 3106 3107 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3108 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3109 PetscFunctionReturn(0); 3110 } 3111 3112 #undef __FUNCT__ 3113 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3114 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3115 { 3116 PetscErrorCode ierr; 3117 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3118 3119 PetscFunctionBegin; 3120 if (coloring->ctype == IS_COLORING_LOCAL) { 3121 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3122 a->coloring = coloring; 3123 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3124 PetscInt i,*larray; 3125 ISColoring ocoloring; 3126 ISColoringValue *colors; 3127 3128 /* set coloring for diagonal portion */ 3129 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3130 for (i=0; i<A->n; i++) { 3131 larray[i] = i; 3132 } 3133 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3134 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3135 for (i=0; i<A->n; i++) { 3136 colors[i] = coloring->colors[larray[i]]; 3137 } 3138 ierr = PetscFree(larray);CHKERRQ(ierr); 3139 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3140 a->coloring = ocoloring; 3141 } 3142 PetscFunctionReturn(0); 3143 } 3144 3145 #if defined(PETSC_HAVE_ADIC) 3146 EXTERN_C_BEGIN 3147 #include "adic/ad_utils.h" 3148 EXTERN_C_END 3149 3150 #undef __FUNCT__ 3151 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3152 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3153 { 3154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3155 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3156 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3157 ISColoringValue *color; 3158 3159 PetscFunctionBegin; 3160 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3161 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3162 color = a->coloring->colors; 3163 /* loop over rows */ 3164 for (i=0; i<m; i++) { 3165 nz = ii[i+1] - ii[i]; 3166 /* loop over columns putting computed value into matrix */ 3167 for (j=0; j<nz; j++) { 3168 *v++ = values[color[*jj++]]; 3169 } 3170 values += nlen; /* jump to next row of derivatives */ 3171 } 3172 PetscFunctionReturn(0); 3173 } 3174 #endif 3175 3176 #undef __FUNCT__ 3177 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3178 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3179 { 3180 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3181 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3182 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3183 ISColoringValue *color; 3184 3185 PetscFunctionBegin; 3186 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3187 color = a->coloring->colors; 3188 /* loop over rows */ 3189 for (i=0; i<m; i++) { 3190 nz = ii[i+1] - ii[i]; 3191 /* loop over columns putting computed value into matrix */ 3192 for (j=0; j<nz; j++) { 3193 *v++ = values[color[*jj++]]; 3194 } 3195 values += nl; /* jump to next row of derivatives */ 3196 } 3197 PetscFunctionReturn(0); 3198 } 3199 3200