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