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 *y; 888 const PetscScalar *x,*aa; 889 PetscErrorCode ierr; 890 PetscInt m=A->rmap.n,*aj,*ii; 891 PetscInt n,i,j,nonzerorow=0,*ridx=PETSC_NULL; 892 PetscScalar sum; 893 PetscTruth usecprow=a->compressedrow.use; 894 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 895 PetscInt jrow; 896 #endif 897 898 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 899 #pragma disjoint(*x,*y,*aa) 900 #endif 901 902 PetscFunctionBegin; 903 ierr = VecGetArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 904 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 905 aj = a->j; 906 aa = a->a; 907 ii = a->i; 908 if (usecprow){ /* use compressed row format */ 909 m = a->compressedrow.nrows; 910 ii = a->compressedrow.i; 911 ridx = a->compressedrow.rindex; 912 for (i=0; i<m; i++){ 913 n = ii[i+1] - ii[i]; 914 aj = a->j + ii[i]; 915 aa = a->a + ii[i]; 916 sum = 0.0; 917 nonzerorow += (n>0); 918 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 919 y[*ridx++] = sum; 920 } 921 } else { /* do not use compressed row format */ 922 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 923 fortranmultaij_(&m,x,ii,aj,aa,y); 924 #else 925 for (i=0; i<m; i++) { 926 jrow = ii[i]; 927 n = ii[i+1] - jrow; 928 sum = 0.0; 929 nonzerorow += (n>0); 930 for (j=0; j<n; j++) { 931 sum += aa[jrow]*x[aj[jrow]]; jrow++; 932 } 933 y[i] = sum; 934 } 935 #endif 936 } 937 ierr = PetscLogFlops(2*a->nz - nonzerorow);CHKERRQ(ierr); 938 ierr = VecRestoreArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 939 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 940 PetscFunctionReturn(0); 941 } 942 943 #undef __FUNCT__ 944 #define __FUNCT__ "MatMultAdd_SeqAIJ" 945 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 946 { 947 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 948 PetscScalar *x,*y,*z,*aa; 949 PetscErrorCode ierr; 950 PetscInt m = A->rmap.n,*aj,*ii; 951 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 952 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 953 PetscScalar sum; 954 PetscTruth usecprow=a->compressedrow.use; 955 #endif 956 957 PetscFunctionBegin; 958 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 959 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 960 if (zz != yy) { 961 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 962 } else { 963 z = y; 964 } 965 966 aj = a->j; 967 aa = a->a; 968 ii = a->i; 969 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 970 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 971 #else 972 if (usecprow){ /* use compressed row format */ 973 if (zz != yy){ 974 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 975 } 976 m = a->compressedrow.nrows; 977 ii = a->compressedrow.i; 978 ridx = a->compressedrow.rindex; 979 for (i=0; i<m; i++){ 980 n = ii[i+1] - ii[i]; 981 aj = a->j + ii[i]; 982 aa = a->a + ii[i]; 983 sum = y[*ridx]; 984 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 985 z[*ridx++] = sum; 986 } 987 } else { /* do not use compressed row format */ 988 for (i=0; i<m; i++) { 989 jrow = ii[i]; 990 n = ii[i+1] - jrow; 991 sum = y[i]; 992 for (j=0; j<n; j++) { 993 sum += aa[jrow]*x[aj[jrow]]; jrow++; 994 } 995 z[i] = sum; 996 } 997 } 998 #endif 999 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1000 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1001 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1002 if (zz != yy) { 1003 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1004 } 1005 PetscFunctionReturn(0); 1006 } 1007 1008 /* 1009 Adds diagonal pointers to sparse matrix structure. 1010 */ 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1013 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1014 { 1015 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1016 PetscErrorCode ierr; 1017 PetscInt i,j,m = A->rmap.n; 1018 1019 PetscFunctionBegin; 1020 if (!a->diag) { 1021 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1022 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1023 } 1024 for (i=0; i<A->rmap.n; i++) { 1025 a->diag[i] = a->i[i+1]; 1026 for (j=a->i[i]; j<a->i[i+1]; j++) { 1027 if (a->j[j] == i) { 1028 a->diag[i] = j; 1029 break; 1030 } 1031 } 1032 } 1033 PetscFunctionReturn(0); 1034 } 1035 1036 /* 1037 Checks for missing diagonals 1038 */ 1039 #undef __FUNCT__ 1040 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1041 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1042 { 1043 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1044 PetscInt *diag,*jj = a->j,i; 1045 1046 PetscFunctionBegin; 1047 *missing = PETSC_FALSE; 1048 if (A->rmap.n > 0 && !jj) { 1049 *missing = PETSC_TRUE; 1050 if (d) *d = 0; 1051 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1052 } else { 1053 diag = a->diag; 1054 for (i=0; i<A->rmap.n; i++) { 1055 if (jj[diag[i]] != i) { 1056 *missing = PETSC_TRUE; 1057 if (d) *d = i; 1058 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1059 } 1060 } 1061 } 1062 PetscFunctionReturn(0); 1063 } 1064 1065 EXTERN_C_BEGIN 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1068 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1069 { 1070 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1071 PetscErrorCode ierr; 1072 PetscInt i,*diag,m = A->rmap.n; 1073 PetscScalar *v = a->a,*idiag,*mdiag; 1074 1075 PetscFunctionBegin; 1076 if (a->idiagvalid) PetscFunctionReturn(0); 1077 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1078 diag = a->diag; 1079 if (!a->idiag) { 1080 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1081 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1082 v = a->a; 1083 } 1084 mdiag = a->mdiag; 1085 idiag = a->idiag; 1086 1087 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1088 for (i=0; i<m; i++) { 1089 mdiag[i] = v[diag[i]]; 1090 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1091 idiag[i] = 1.0/v[diag[i]]; 1092 } 1093 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1094 } else { 1095 for (i=0; i<m; i++) { 1096 mdiag[i] = v[diag[i]]; 1097 idiag[i] = omega/(fshift + v[diag[i]]); 1098 } 1099 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1100 } 1101 a->idiagvalid = PETSC_TRUE; 1102 PetscFunctionReturn(0); 1103 } 1104 EXTERN_C_END 1105 1106 #undef __FUNCT__ 1107 #define __FUNCT__ "MatRelax_SeqAIJ" 1108 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1109 { 1110 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1111 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1112 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1113 PetscErrorCode ierr; 1114 PetscInt n = A->cmap.n,m = A->rmap.n,i; 1115 const PetscInt *idx,*diag; 1116 1117 PetscFunctionBegin; 1118 its = its*lits; 1119 1120 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1121 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1122 a->fshift = fshift; 1123 a->omega = omega; 1124 1125 diag = a->diag; 1126 t = a->ssor_work; 1127 idiag = a->idiag; 1128 mdiag = a->mdiag; 1129 1130 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1131 if (xx != bb) { 1132 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1133 } else { 1134 b = x; 1135 } 1136 CHKMEMQ; 1137 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1138 xs = x; 1139 if (flag == SOR_APPLY_UPPER) { 1140 /* apply (U + D/omega) to the vector */ 1141 bs = b; 1142 for (i=0; i<m; i++) { 1143 d = fshift + mdiag[i]; 1144 n = a->i[i+1] - diag[i] - 1; 1145 idx = a->j + diag[i] + 1; 1146 v = a->a + diag[i] + 1; 1147 sum = b[i]*d/omega; 1148 SPARSEDENSEDOT(sum,bs,v,idx,n); 1149 x[i] = sum; 1150 } 1151 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1152 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1153 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1154 PetscFunctionReturn(0); 1155 } 1156 1157 if (flag == SOR_APPLY_LOWER) { 1158 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1159 } else if (flag & SOR_EISENSTAT) { 1160 /* Let A = L + U + D; where L is lower trianglar, 1161 U is upper triangular, E is diagonal; This routine applies 1162 1163 (L + E)^{-1} A (U + E)^{-1} 1164 1165 to a vector efficiently using Eisenstat's trick. This is for 1166 the case of SSOR preconditioner, so E is D/omega where omega 1167 is the relaxation factor. 1168 */ 1169 scale = (2.0/omega) - 1.0; 1170 1171 /* x = (E + U)^{-1} b */ 1172 for (i=m-1; i>=0; i--) { 1173 n = a->i[i+1] - diag[i] - 1; 1174 idx = a->j + diag[i] + 1; 1175 v = a->a + diag[i] + 1; 1176 sum = b[i]; 1177 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1178 x[i] = sum*idiag[i]; 1179 } 1180 1181 /* t = b - (2*E - D)x */ 1182 v = a->a; 1183 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1184 1185 /* t = (E + L)^{-1}t */ 1186 ts = t; 1187 diag = a->diag; 1188 for (i=0; i<m; i++) { 1189 n = diag[i] - a->i[i]; 1190 idx = a->j + a->i[i]; 1191 v = a->a + a->i[i]; 1192 sum = t[i]; 1193 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1194 t[i] = sum*idiag[i]; 1195 /* x = x + t */ 1196 x[i] += t[i]; 1197 } 1198 1199 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1200 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1201 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1202 PetscFunctionReturn(0); 1203 } 1204 if (flag & SOR_ZERO_INITIAL_GUESS) { 1205 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1206 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1207 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1208 #else 1209 for (i=0; i<m; i++) { 1210 n = diag[i] - a->i[i]; 1211 idx = a->j + a->i[i]; 1212 v = a->a + a->i[i]; 1213 sum = b[i]; 1214 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1215 x[i] = sum*idiag[i]; 1216 } 1217 #endif 1218 xb = x; 1219 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1220 } else xb = b; 1221 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1222 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1223 for (i=0; i<m; i++) { 1224 x[i] *= mdiag[i]; 1225 } 1226 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1227 } 1228 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1229 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1230 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1231 #else 1232 for (i=m-1; i>=0; i--) { 1233 n = a->i[i+1] - diag[i] - 1; 1234 idx = a->j + diag[i] + 1; 1235 v = a->a + diag[i] + 1; 1236 sum = xb[i]; 1237 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1238 x[i] = sum*idiag[i]; 1239 } 1240 #endif 1241 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1242 } 1243 its--; 1244 } 1245 while (its--) { 1246 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1247 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1248 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1249 #else 1250 for (i=0; i<m; i++) { 1251 n = a->i[i+1] - a->i[i]; 1252 idx = a->j + a->i[i]; 1253 v = a->a + a->i[i]; 1254 sum = b[i]; 1255 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1256 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1257 } 1258 #endif 1259 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1260 } 1261 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1262 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1263 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1264 #else 1265 for (i=m-1; i>=0; i--) { 1266 n = a->i[i+1] - a->i[i]; 1267 idx = a->j + a->i[i]; 1268 v = a->a + a->i[i]; 1269 sum = b[i]; 1270 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1271 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1272 } 1273 #endif 1274 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1275 } 1276 } 1277 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1278 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1279 CHKMEMQ; PetscFunctionReturn(0); 1280 } 1281 1282 1283 #undef __FUNCT__ 1284 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1285 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1286 { 1287 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1288 1289 PetscFunctionBegin; 1290 info->rows_global = (double)A->rmap.n; 1291 info->columns_global = (double)A->cmap.n; 1292 info->rows_local = (double)A->rmap.n; 1293 info->columns_local = (double)A->cmap.n; 1294 info->block_size = 1.0; 1295 info->nz_allocated = (double)a->maxnz; 1296 info->nz_used = (double)a->nz; 1297 info->nz_unneeded = (double)(a->maxnz - a->nz); 1298 info->assemblies = (double)A->num_ass; 1299 info->mallocs = (double)a->reallocs; 1300 info->memory = ((PetscObject)A)->mem; 1301 if (A->factor) { 1302 info->fill_ratio_given = A->info.fill_ratio_given; 1303 info->fill_ratio_needed = A->info.fill_ratio_needed; 1304 info->factor_mallocs = A->info.factor_mallocs; 1305 } else { 1306 info->fill_ratio_given = 0; 1307 info->fill_ratio_needed = 0; 1308 info->factor_mallocs = 0; 1309 } 1310 PetscFunctionReturn(0); 1311 } 1312 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1315 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1316 { 1317 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1318 PetscInt i,m = A->rmap.n - 1,d; 1319 PetscErrorCode ierr; 1320 PetscTruth missing; 1321 1322 PetscFunctionBegin; 1323 if (a->keepzeroedrows) { 1324 for (i=0; i<N; i++) { 1325 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1326 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1327 } 1328 if (diag != 0.0) { 1329 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1330 if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1331 for (i=0; i<N; i++) { 1332 a->a[a->diag[rows[i]]] = diag; 1333 } 1334 } 1335 A->same_nonzero = PETSC_TRUE; 1336 } else { 1337 if (diag != 0.0) { 1338 for (i=0; i<N; i++) { 1339 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1340 if (a->ilen[rows[i]] > 0) { 1341 a->ilen[rows[i]] = 1; 1342 a->a[a->i[rows[i]]] = diag; 1343 a->j[a->i[rows[i]]] = rows[i]; 1344 } else { /* in case row was completely empty */ 1345 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1346 } 1347 } 1348 } else { 1349 for (i=0; i<N; i++) { 1350 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1351 a->ilen[rows[i]] = 0; 1352 } 1353 } 1354 A->same_nonzero = PETSC_FALSE; 1355 } 1356 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1357 PetscFunctionReturn(0); 1358 } 1359 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "MatGetRow_SeqAIJ" 1362 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1363 { 1364 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1365 PetscInt *itmp; 1366 1367 PetscFunctionBegin; 1368 if (row < 0 || row >= A->rmap.n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1369 1370 *nz = a->i[row+1] - a->i[row]; 1371 if (v) *v = a->a + a->i[row]; 1372 if (idx) { 1373 itmp = a->j + a->i[row]; 1374 if (*nz) { 1375 *idx = itmp; 1376 } 1377 else *idx = 0; 1378 } 1379 PetscFunctionReturn(0); 1380 } 1381 1382 /* remove this function? */ 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1385 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1386 { 1387 PetscFunctionBegin; 1388 PetscFunctionReturn(0); 1389 } 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "MatNorm_SeqAIJ" 1393 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1394 { 1395 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1396 PetscScalar *v = a->a; 1397 PetscReal sum = 0.0; 1398 PetscErrorCode ierr; 1399 PetscInt i,j; 1400 1401 PetscFunctionBegin; 1402 if (type == NORM_FROBENIUS) { 1403 for (i=0; i<a->nz; i++) { 1404 #if defined(PETSC_USE_COMPLEX) 1405 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1406 #else 1407 sum += (*v)*(*v); v++; 1408 #endif 1409 } 1410 *nrm = sqrt(sum); 1411 } else if (type == NORM_1) { 1412 PetscReal *tmp; 1413 PetscInt *jj = a->j; 1414 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1415 ierr = PetscMemzero(tmp,A->cmap.n*sizeof(PetscReal));CHKERRQ(ierr); 1416 *nrm = 0.0; 1417 for (j=0; j<a->nz; j++) { 1418 tmp[*jj++] += PetscAbsScalar(*v); v++; 1419 } 1420 for (j=0; j<A->cmap.n; j++) { 1421 if (tmp[j] > *nrm) *nrm = tmp[j]; 1422 } 1423 ierr = PetscFree(tmp);CHKERRQ(ierr); 1424 } else if (type == NORM_INFINITY) { 1425 *nrm = 0.0; 1426 for (j=0; j<A->rmap.n; j++) { 1427 v = a->a + a->i[j]; 1428 sum = 0.0; 1429 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1430 sum += PetscAbsScalar(*v); v++; 1431 } 1432 if (sum > *nrm) *nrm = sum; 1433 } 1434 } else { 1435 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1436 } 1437 PetscFunctionReturn(0); 1438 } 1439 1440 #undef __FUNCT__ 1441 #define __FUNCT__ "MatTranspose_SeqAIJ" 1442 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1443 { 1444 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1445 Mat C; 1446 PetscErrorCode ierr; 1447 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap.n,len,*col; 1448 PetscScalar *array = a->a; 1449 1450 PetscFunctionBegin; 1451 if (!B && m != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1452 ierr = PetscMalloc((1+A->cmap.n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1453 ierr = PetscMemzero(col,(1+A->cmap.n)*sizeof(PetscInt));CHKERRQ(ierr); 1454 1455 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1456 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1457 ierr = MatSetSizes(C,A->cmap.n,m,A->cmap.n,m);CHKERRQ(ierr); 1458 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1459 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1460 ierr = PetscFree(col);CHKERRQ(ierr); 1461 for (i=0; i<m; i++) { 1462 len = ai[i+1]-ai[i]; 1463 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1464 array += len; 1465 aj += len; 1466 } 1467 1468 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1469 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1470 1471 if (B) { 1472 *B = C; 1473 } else { 1474 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1475 } 1476 PetscFunctionReturn(0); 1477 } 1478 1479 EXTERN_C_BEGIN 1480 #undef __FUNCT__ 1481 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1482 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1483 { 1484 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1485 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1486 PetscErrorCode ierr; 1487 PetscInt ma,na,mb,nb, i; 1488 1489 PetscFunctionBegin; 1490 bij = (Mat_SeqAIJ *) B->data; 1491 1492 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1493 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1494 if (ma!=nb || na!=mb){ 1495 *f = PETSC_FALSE; 1496 PetscFunctionReturn(0); 1497 } 1498 aii = aij->i; bii = bij->i; 1499 adx = aij->j; bdx = bij->j; 1500 va = aij->a; vb = bij->a; 1501 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1502 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1503 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1504 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1505 1506 *f = PETSC_TRUE; 1507 for (i=0; i<ma; i++) { 1508 while (aptr[i]<aii[i+1]) { 1509 PetscInt idc,idr; 1510 PetscScalar vc,vr; 1511 /* column/row index/value */ 1512 idc = adx[aptr[i]]; 1513 idr = bdx[bptr[idc]]; 1514 vc = va[aptr[i]]; 1515 vr = vb[bptr[idc]]; 1516 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1517 *f = PETSC_FALSE; 1518 goto done; 1519 } else { 1520 aptr[i]++; 1521 if (B || i!=idc) bptr[idc]++; 1522 } 1523 } 1524 } 1525 done: 1526 ierr = PetscFree(aptr);CHKERRQ(ierr); 1527 if (B) { 1528 ierr = PetscFree(bptr);CHKERRQ(ierr); 1529 } 1530 PetscFunctionReturn(0); 1531 } 1532 EXTERN_C_END 1533 1534 EXTERN_C_BEGIN 1535 #undef __FUNCT__ 1536 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1537 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1538 { 1539 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1540 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1541 PetscErrorCode ierr; 1542 PetscInt ma,na,mb,nb, i; 1543 1544 PetscFunctionBegin; 1545 bij = (Mat_SeqAIJ *) B->data; 1546 1547 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1548 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1549 if (ma!=nb || na!=mb){ 1550 *f = PETSC_FALSE; 1551 PetscFunctionReturn(0); 1552 } 1553 aii = aij->i; bii = bij->i; 1554 adx = aij->j; bdx = bij->j; 1555 va = aij->a; vb = bij->a; 1556 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1557 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1558 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1559 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1560 1561 *f = PETSC_TRUE; 1562 for (i=0; i<ma; i++) { 1563 while (aptr[i]<aii[i+1]) { 1564 PetscInt idc,idr; 1565 PetscScalar vc,vr; 1566 /* column/row index/value */ 1567 idc = adx[aptr[i]]; 1568 idr = bdx[bptr[idc]]; 1569 vc = va[aptr[i]]; 1570 vr = vb[bptr[idc]]; 1571 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1572 *f = PETSC_FALSE; 1573 goto done; 1574 } else { 1575 aptr[i]++; 1576 if (B || i!=idc) bptr[idc]++; 1577 } 1578 } 1579 } 1580 done: 1581 ierr = PetscFree(aptr);CHKERRQ(ierr); 1582 if (B) { 1583 ierr = PetscFree(bptr);CHKERRQ(ierr); 1584 } 1585 PetscFunctionReturn(0); 1586 } 1587 EXTERN_C_END 1588 1589 #undef __FUNCT__ 1590 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1591 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1592 { 1593 PetscErrorCode ierr; 1594 PetscFunctionBegin; 1595 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1601 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1602 { 1603 PetscErrorCode ierr; 1604 PetscFunctionBegin; 1605 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1611 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1612 { 1613 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1614 PetscScalar *l,*r,x,*v; 1615 PetscErrorCode ierr; 1616 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,M,nz = a->nz,*jj; 1617 1618 PetscFunctionBegin; 1619 if (ll) { 1620 /* The local size is used so that VecMPI can be passed to this routine 1621 by MatDiagonalScale_MPIAIJ */ 1622 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1623 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1624 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1625 v = a->a; 1626 for (i=0; i<m; i++) { 1627 x = l[i]; 1628 M = a->i[i+1] - a->i[i]; 1629 for (j=0; j<M; j++) { (*v++) *= x;} 1630 } 1631 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1632 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1633 } 1634 if (rr) { 1635 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1636 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1637 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1638 v = a->a; jj = a->j; 1639 for (i=0; i<nz; i++) { 1640 (*v++) *= r[*jj++]; 1641 } 1642 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1643 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1644 } 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #undef __FUNCT__ 1649 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1650 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1651 { 1652 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1653 PetscErrorCode ierr; 1654 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap.n,*lens; 1655 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1656 PetscInt *irow,*icol,nrows,ncols; 1657 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1658 PetscScalar *a_new,*mat_a; 1659 Mat C; 1660 PetscTruth stride; 1661 1662 PetscFunctionBegin; 1663 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1664 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1665 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1666 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1667 1668 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1669 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1670 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1671 1672 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1673 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1674 if (stride && step == 1) { 1675 /* special case of contiguous rows */ 1676 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1677 starts = lens + nrows; 1678 /* loop over new rows determining lens and starting points */ 1679 for (i=0; i<nrows; i++) { 1680 kstart = ai[irow[i]]; 1681 kend = kstart + ailen[irow[i]]; 1682 for (k=kstart; k<kend; k++) { 1683 if (aj[k] >= first) { 1684 starts[i] = k; 1685 break; 1686 } 1687 } 1688 sum = 0; 1689 while (k < kend) { 1690 if (aj[k++] >= first+ncols) break; 1691 sum++; 1692 } 1693 lens[i] = sum; 1694 } 1695 /* create submatrix */ 1696 if (scall == MAT_REUSE_MATRIX) { 1697 PetscInt n_cols,n_rows; 1698 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1699 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1700 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1701 C = *B; 1702 } else { 1703 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1704 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1705 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1706 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1707 } 1708 c = (Mat_SeqAIJ*)C->data; 1709 1710 /* loop over rows inserting into submatrix */ 1711 a_new = c->a; 1712 j_new = c->j; 1713 i_new = c->i; 1714 1715 for (i=0; i<nrows; i++) { 1716 ii = starts[i]; 1717 lensi = lens[i]; 1718 for (k=0; k<lensi; k++) { 1719 *j_new++ = aj[ii+k] - first; 1720 } 1721 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1722 a_new += lensi; 1723 i_new[i+1] = i_new[i] + lensi; 1724 c->ilen[i] = lensi; 1725 } 1726 ierr = PetscFree(lens);CHKERRQ(ierr); 1727 } else { 1728 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1729 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1730 1731 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1732 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1733 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1734 /* determine lens of each row */ 1735 for (i=0; i<nrows; i++) { 1736 kstart = ai[irow[i]]; 1737 kend = kstart + a->ilen[irow[i]]; 1738 lens[i] = 0; 1739 for (k=kstart; k<kend; k++) { 1740 if (smap[aj[k]]) { 1741 lens[i]++; 1742 } 1743 } 1744 } 1745 /* Create and fill new matrix */ 1746 if (scall == MAT_REUSE_MATRIX) { 1747 PetscTruth equal; 1748 1749 c = (Mat_SeqAIJ *)((*B)->data); 1750 if ((*B)->rmap.n != nrows || (*B)->cmap.n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1751 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap.n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1752 if (!equal) { 1753 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1754 } 1755 ierr = PetscMemzero(c->ilen,(*B)->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 1756 C = *B; 1757 } else { 1758 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1759 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1760 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1761 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1762 } 1763 c = (Mat_SeqAIJ *)(C->data); 1764 for (i=0; i<nrows; i++) { 1765 row = irow[i]; 1766 kstart = ai[row]; 1767 kend = kstart + a->ilen[row]; 1768 mat_i = c->i[i]; 1769 mat_j = c->j + mat_i; 1770 mat_a = c->a + mat_i; 1771 mat_ilen = c->ilen + i; 1772 for (k=kstart; k<kend; k++) { 1773 if ((tcol=smap[a->j[k]])) { 1774 *mat_j++ = tcol - 1; 1775 *mat_a++ = a->a[k]; 1776 (*mat_ilen)++; 1777 1778 } 1779 } 1780 } 1781 /* Free work space */ 1782 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1783 ierr = PetscFree(smap);CHKERRQ(ierr); 1784 ierr = PetscFree(lens);CHKERRQ(ierr); 1785 } 1786 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1787 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1788 1789 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1790 *B = C; 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /* 1795 */ 1796 #undef __FUNCT__ 1797 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1798 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1799 { 1800 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1801 PetscErrorCode ierr; 1802 Mat outA; 1803 PetscTruth row_identity,col_identity; 1804 1805 PetscFunctionBegin; 1806 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1807 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1808 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1809 1810 outA = inA; 1811 inA->factor = FACTOR_LU; 1812 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1813 if (a->row) { ierr = ISDestroy(a->row); CHKERRQ(ierr);} 1814 a->row = row; 1815 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1816 if (a->col) { ierr = ISDestroy(a->col); CHKERRQ(ierr);} 1817 a->col = col; 1818 1819 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1820 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1821 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1822 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1823 1824 if (!a->solve_work) { /* this matrix may have been factored before */ 1825 ierr = PetscMalloc((inA->rmap.n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1826 ierr = PetscLogObjectMemory(inA, (inA->rmap.n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1827 } 1828 1829 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1830 if (row_identity && col_identity) { 1831 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1832 } else { 1833 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(inA,info,&outA);CHKERRQ(ierr); 1834 } 1835 PetscFunctionReturn(0); 1836 } 1837 1838 #include "petscblaslapack.h" 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "MatScale_SeqAIJ" 1841 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1842 { 1843 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1844 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1845 PetscScalar oalpha = alpha; 1846 PetscErrorCode ierr; 1847 1848 1849 PetscFunctionBegin; 1850 BLASscal_(&bnz,&oalpha,a->a,&one); 1851 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1857 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1858 { 1859 PetscErrorCode ierr; 1860 PetscInt i; 1861 1862 PetscFunctionBegin; 1863 if (scall == MAT_INITIAL_MATRIX) { 1864 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1865 } 1866 1867 for (i=0; i<n; i++) { 1868 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 #undef __FUNCT__ 1874 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1875 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1876 { 1877 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1878 PetscErrorCode ierr; 1879 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1880 PetscInt start,end,*ai,*aj; 1881 PetscBT table; 1882 1883 PetscFunctionBegin; 1884 m = A->rmap.n; 1885 ai = a->i; 1886 aj = a->j; 1887 1888 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1889 1890 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1891 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1892 1893 for (i=0; i<is_max; i++) { 1894 /* Initialize the two local arrays */ 1895 isz = 0; 1896 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1897 1898 /* Extract the indices, assume there can be duplicate entries */ 1899 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1900 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1901 1902 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1903 for (j=0; j<n ; ++j){ 1904 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1905 } 1906 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1907 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1908 1909 k = 0; 1910 for (j=0; j<ov; j++){ /* for each overlap */ 1911 n = isz; 1912 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1913 row = nidx[k]; 1914 start = ai[row]; 1915 end = ai[row+1]; 1916 for (l = start; l<end ; l++){ 1917 val = aj[l] ; 1918 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1919 } 1920 } 1921 } 1922 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1923 } 1924 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1925 ierr = PetscFree(nidx);CHKERRQ(ierr); 1926 PetscFunctionReturn(0); 1927 } 1928 1929 /* -------------------------------------------------------------- */ 1930 #undef __FUNCT__ 1931 #define __FUNCT__ "MatPermute_SeqAIJ" 1932 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1933 { 1934 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1935 PetscErrorCode ierr; 1936 PetscInt i,nz,m = A->rmap.n,n = A->cmap.n,*col; 1937 PetscInt *row,*cnew,j,*lens; 1938 IS icolp,irowp; 1939 PetscInt *cwork; 1940 PetscScalar *vwork; 1941 1942 PetscFunctionBegin; 1943 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1944 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1945 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1946 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1947 1948 /* determine lengths of permuted rows */ 1949 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1950 for (i=0; i<m; i++) { 1951 lens[row[i]] = a->i[i+1] - a->i[i]; 1952 } 1953 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1954 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1955 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1956 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1957 ierr = PetscFree(lens);CHKERRQ(ierr); 1958 1959 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1960 for (i=0; i<m; i++) { 1961 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1962 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1963 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1964 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1965 } 1966 ierr = PetscFree(cnew);CHKERRQ(ierr); 1967 (*B)->assembled = PETSC_FALSE; 1968 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1969 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1970 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1971 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1972 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1973 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1974 PetscFunctionReturn(0); 1975 } 1976 1977 #undef __FUNCT__ 1978 #define __FUNCT__ "MatCopy_SeqAIJ" 1979 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1980 { 1981 PetscErrorCode ierr; 1982 1983 PetscFunctionBegin; 1984 /* If the two matrices have the same copy implementation, use fast copy. */ 1985 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1986 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1987 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1988 1989 if (a->i[A->rmap.n] != b->i[B->rmap.n]) { 1990 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1991 } 1992 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 1993 } else { 1994 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1995 } 1996 PetscFunctionReturn(0); 1997 } 1998 1999 #undef __FUNCT__ 2000 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2001 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2002 { 2003 PetscErrorCode ierr; 2004 2005 PetscFunctionBegin; 2006 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2007 PetscFunctionReturn(0); 2008 } 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "MatGetArray_SeqAIJ" 2012 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2013 { 2014 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2015 PetscFunctionBegin; 2016 *array = a->a; 2017 PetscFunctionReturn(0); 2018 } 2019 2020 #undef __FUNCT__ 2021 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2022 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2023 { 2024 PetscFunctionBegin; 2025 PetscFunctionReturn(0); 2026 } 2027 2028 #undef __FUNCT__ 2029 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2030 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2031 { 2032 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2033 PetscErrorCode ierr; 2034 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2035 PetscScalar dx,*y,*xx,*w3_array; 2036 PetscScalar *vscale_array; 2037 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2038 Vec w1,w2,w3; 2039 void *fctx = coloring->fctx; 2040 PetscTruth flg; 2041 2042 PetscFunctionBegin; 2043 if (!coloring->w1) { 2044 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2045 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2046 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2047 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2048 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2049 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2050 } 2051 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2052 2053 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2054 ierr = PetscOptionsHasName(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 2055 if (flg) { 2056 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2057 } else { 2058 PetscTruth assembled; 2059 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2060 if (assembled) { 2061 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2062 } 2063 } 2064 2065 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2066 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2067 2068 /* 2069 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2070 coloring->F for the coarser grids from the finest 2071 */ 2072 if (coloring->F) { 2073 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2074 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2075 if (m1 != m2) { 2076 coloring->F = 0; 2077 } 2078 } 2079 2080 if (coloring->F) { 2081 w1 = coloring->F; 2082 coloring->F = 0; 2083 } else { 2084 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2085 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2086 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2087 } 2088 2089 /* 2090 Compute all the scale factors and share with other processors 2091 */ 2092 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2093 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2094 for (k=0; k<coloring->ncolors; k++) { 2095 /* 2096 Loop over each column associated with color adding the 2097 perturbation to the vector w3. 2098 */ 2099 for (l=0; l<coloring->ncolumns[k]; l++) { 2100 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2101 dx = xx[col]; 2102 if (dx == 0.0) dx = 1.0; 2103 #if !defined(PETSC_USE_COMPLEX) 2104 if (dx < umin && dx >= 0.0) dx = umin; 2105 else if (dx < 0.0 && dx > -umin) dx = -umin; 2106 #else 2107 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2108 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2109 #endif 2110 dx *= epsilon; 2111 vscale_array[col] = 1.0/dx; 2112 } 2113 } 2114 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2115 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2116 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2117 2118 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2119 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2120 2121 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2122 else vscaleforrow = coloring->columnsforrow; 2123 2124 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2125 /* 2126 Loop over each color 2127 */ 2128 for (k=0; k<coloring->ncolors; k++) { 2129 coloring->currentcolor = k; 2130 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2131 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2132 /* 2133 Loop over each column associated with color adding the 2134 perturbation to the vector w3. 2135 */ 2136 for (l=0; l<coloring->ncolumns[k]; l++) { 2137 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2138 dx = xx[col]; 2139 if (dx == 0.0) dx = 1.0; 2140 #if !defined(PETSC_USE_COMPLEX) 2141 if (dx < umin && dx >= 0.0) dx = umin; 2142 else if (dx < 0.0 && dx > -umin) dx = -umin; 2143 #else 2144 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2145 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2146 #endif 2147 dx *= epsilon; 2148 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2149 w3_array[col] += dx; 2150 } 2151 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2152 2153 /* 2154 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2155 */ 2156 2157 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2158 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2159 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2160 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2161 2162 /* 2163 Loop over rows of vector, putting results into Jacobian matrix 2164 */ 2165 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2166 for (l=0; l<coloring->nrows[k]; l++) { 2167 row = coloring->rows[k][l]; 2168 col = coloring->columnsforrow[k][l]; 2169 y[row] *= vscale_array[vscaleforrow[k][l]]; 2170 srow = row + start; 2171 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2172 } 2173 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2174 } 2175 coloring->currentcolor = k; 2176 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2177 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2178 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2179 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2180 PetscFunctionReturn(0); 2181 } 2182 2183 #include "petscblaslapack.h" 2184 #undef __FUNCT__ 2185 #define __FUNCT__ "MatAXPY_SeqAIJ" 2186 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2187 { 2188 PetscErrorCode ierr; 2189 PetscInt i; 2190 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2191 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2192 2193 PetscFunctionBegin; 2194 if (str == SAME_NONZERO_PATTERN) { 2195 PetscScalar alpha = a; 2196 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2197 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2198 if (y->xtoy && y->XtoY != X) { 2199 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2200 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2201 } 2202 if (!y->xtoy) { /* get xtoy */ 2203 ierr = MatAXPYGetxtoy_Private(X->rmap.n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2204 y->XtoY = X; 2205 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2206 } 2207 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2208 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2209 } else { 2210 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2211 } 2212 PetscFunctionReturn(0); 2213 } 2214 2215 #undef __FUNCT__ 2216 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2217 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2218 { 2219 PetscFunctionBegin; 2220 PetscFunctionReturn(0); 2221 } 2222 2223 #undef __FUNCT__ 2224 #define __FUNCT__ "MatConjugate_SeqAIJ" 2225 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2226 { 2227 #if defined(PETSC_USE_COMPLEX) 2228 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2229 PetscInt i,nz; 2230 PetscScalar *a; 2231 2232 PetscFunctionBegin; 2233 nz = aij->nz; 2234 a = aij->a; 2235 for (i=0; i<nz; i++) { 2236 a[i] = PetscConj(a[i]); 2237 } 2238 #else 2239 PetscFunctionBegin; 2240 #endif 2241 PetscFunctionReturn(0); 2242 } 2243 2244 #undef __FUNCT__ 2245 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2246 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2247 { 2248 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2249 PetscErrorCode ierr; 2250 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2251 PetscReal atmp; 2252 PetscScalar *x; 2253 MatScalar *aa; 2254 2255 PetscFunctionBegin; 2256 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2257 aa = a->a; 2258 ai = a->i; 2259 aj = a->j; 2260 2261 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2262 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2263 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2264 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2265 for (i=0; i<m; i++) { 2266 ncols = ai[1] - ai[0]; ai++; 2267 x[i] = 0.0; if (idx) idx[i] = 0; 2268 for (j=0; j<ncols; j++){ 2269 atmp = PetscAbsScalar(*aa); 2270 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2271 aa++; aj++; 2272 } 2273 } 2274 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2275 PetscFunctionReturn(0); 2276 } 2277 2278 #undef __FUNCT__ 2279 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2280 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2281 { 2282 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2283 PetscErrorCode ierr; 2284 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2285 PetscScalar *x; 2286 MatScalar *aa; 2287 2288 PetscFunctionBegin; 2289 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2290 aa = a->a; 2291 ai = a->i; 2292 aj = a->j; 2293 2294 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2295 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2296 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2297 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2298 for (i=0; i<m; i++) { 2299 ncols = ai[1] - ai[0]; ai++; 2300 if (ncols == A->cmap.n) { /* row is dense */ 2301 x[i] = *aa; if (idx) idx[i] = 0; 2302 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2303 x[i] = 0.0; 2304 if (idx) { 2305 idx[i] = 0; /* in case ncols is zero */ 2306 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2307 if (aj[j] > j) { 2308 idx[i] = j; 2309 break; 2310 } 2311 } 2312 } 2313 } 2314 for (j=0; j<ncols; j++){ 2315 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2316 aa++; aj++; 2317 } 2318 } 2319 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 #undef __FUNCT__ 2324 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2325 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2326 { 2327 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2328 PetscErrorCode ierr; 2329 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2330 PetscScalar *x; 2331 MatScalar *aa; 2332 2333 PetscFunctionBegin; 2334 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2335 aa = a->a; 2336 ai = a->i; 2337 aj = a->j; 2338 2339 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2340 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2341 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2342 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2343 for (i=0; i<m; i++) { 2344 ncols = ai[1] - ai[0]; ai++; 2345 if (ncols == A->cmap.n) { /* row is dense */ 2346 x[i] = *aa; if (idx) idx[i] = 0; 2347 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2348 x[i] = 0.0; 2349 if (idx) { /* find first implicit 0.0 in the row */ 2350 idx[i] = 0; /* in case ncols is zero */ 2351 for (j=0;j<ncols;j++) { 2352 if (aj[j] > j) { 2353 idx[i] = j; 2354 break; 2355 } 2356 } 2357 } 2358 } 2359 for (j=0; j<ncols; j++){ 2360 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2361 aa++; aj++; 2362 } 2363 } 2364 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 /* -------------------------------------------------------------------*/ 2369 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2370 MatGetRow_SeqAIJ, 2371 MatRestoreRow_SeqAIJ, 2372 MatMult_SeqAIJ, 2373 /* 4*/ MatMultAdd_SeqAIJ, 2374 MatMultTranspose_SeqAIJ, 2375 MatMultTransposeAdd_SeqAIJ, 2376 MatSolve_SeqAIJ, 2377 MatSolveAdd_SeqAIJ, 2378 MatSolveTranspose_SeqAIJ, 2379 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2380 MatLUFactor_SeqAIJ, 2381 0, 2382 MatRelax_SeqAIJ, 2383 MatTranspose_SeqAIJ, 2384 /*15*/ MatGetInfo_SeqAIJ, 2385 MatEqual_SeqAIJ, 2386 MatGetDiagonal_SeqAIJ, 2387 MatDiagonalScale_SeqAIJ, 2388 MatNorm_SeqAIJ, 2389 /*20*/ 0, 2390 MatAssemblyEnd_SeqAIJ, 2391 MatCompress_SeqAIJ, 2392 MatSetOption_SeqAIJ, 2393 MatZeroEntries_SeqAIJ, 2394 /*25*/ MatZeroRows_SeqAIJ, 2395 MatLUFactorSymbolic_SeqAIJ, 2396 MatLUFactorNumeric_SeqAIJ, 2397 MatCholeskyFactorSymbolic_SeqAIJ, 2398 MatCholeskyFactorNumeric_SeqAIJ, 2399 /*30*/ MatSetUpPreallocation_SeqAIJ, 2400 MatILUFactorSymbolic_SeqAIJ, 2401 MatICCFactorSymbolic_SeqAIJ, 2402 MatGetArray_SeqAIJ, 2403 MatRestoreArray_SeqAIJ, 2404 /*35*/ MatDuplicate_SeqAIJ, 2405 0, 2406 0, 2407 MatILUFactor_SeqAIJ, 2408 0, 2409 /*40*/ MatAXPY_SeqAIJ, 2410 MatGetSubMatrices_SeqAIJ, 2411 MatIncreaseOverlap_SeqAIJ, 2412 MatGetValues_SeqAIJ, 2413 MatCopy_SeqAIJ, 2414 /*45*/ MatGetRowMax_SeqAIJ, 2415 MatScale_SeqAIJ, 2416 0, 2417 MatDiagonalSet_SeqAIJ, 2418 MatILUDTFactor_SeqAIJ, 2419 /*50*/ MatSetBlockSize_SeqAIJ, 2420 MatGetRowIJ_SeqAIJ, 2421 MatRestoreRowIJ_SeqAIJ, 2422 MatGetColumnIJ_SeqAIJ, 2423 MatRestoreColumnIJ_SeqAIJ, 2424 /*55*/ MatFDColoringCreate_SeqAIJ, 2425 0, 2426 0, 2427 MatPermute_SeqAIJ, 2428 0, 2429 /*60*/ 0, 2430 MatDestroy_SeqAIJ, 2431 MatView_SeqAIJ, 2432 0, 2433 0, 2434 /*65*/ 0, 2435 0, 2436 0, 2437 0, 2438 0, 2439 /*70*/ MatGetRowMaxAbs_SeqAIJ, 2440 0, 2441 MatSetColoring_SeqAIJ, 2442 #if defined(PETSC_HAVE_ADIC) 2443 MatSetValuesAdic_SeqAIJ, 2444 #else 2445 0, 2446 #endif 2447 MatSetValuesAdifor_SeqAIJ, 2448 /*75*/ 0, 2449 0, 2450 0, 2451 0, 2452 0, 2453 /*80*/ 0, 2454 0, 2455 0, 2456 0, 2457 MatLoad_SeqAIJ, 2458 /*85*/ MatIsSymmetric_SeqAIJ, 2459 MatIsHermitian_SeqAIJ, 2460 0, 2461 0, 2462 0, 2463 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2464 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2465 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2466 MatPtAP_Basic, 2467 MatPtAPSymbolic_SeqAIJ, 2468 /*95*/ MatPtAPNumeric_SeqAIJ, 2469 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2470 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2471 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2472 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2473 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2474 0, 2475 0, 2476 MatConjugate_SeqAIJ, 2477 0, 2478 /*105*/MatSetValuesRow_SeqAIJ, 2479 MatRealPart_SeqAIJ, 2480 MatImaginaryPart_SeqAIJ, 2481 0, 2482 0, 2483 /*110*/MatMatSolve_SeqAIJ, 2484 0, 2485 MatGetRowMin_SeqAIJ, 2486 0, 2487 MatMissingDiagonal_SeqAIJ 2488 }; 2489 2490 EXTERN_C_BEGIN 2491 #undef __FUNCT__ 2492 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2493 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2494 { 2495 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2496 PetscInt i,nz,n; 2497 2498 PetscFunctionBegin; 2499 2500 nz = aij->maxnz; 2501 n = mat->rmap.n; 2502 for (i=0; i<nz; i++) { 2503 aij->j[i] = indices[i]; 2504 } 2505 aij->nz = nz; 2506 for (i=0; i<n; i++) { 2507 aij->ilen[i] = aij->imax[i]; 2508 } 2509 2510 PetscFunctionReturn(0); 2511 } 2512 EXTERN_C_END 2513 2514 #undef __FUNCT__ 2515 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2516 /*@ 2517 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2518 in the matrix. 2519 2520 Input Parameters: 2521 + mat - the SeqAIJ matrix 2522 - indices - the column indices 2523 2524 Level: advanced 2525 2526 Notes: 2527 This can be called if you have precomputed the nonzero structure of the 2528 matrix and want to provide it to the matrix object to improve the performance 2529 of the MatSetValues() operation. 2530 2531 You MUST have set the correct numbers of nonzeros per row in the call to 2532 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2533 2534 MUST be called before any calls to MatSetValues(); 2535 2536 The indices should start with zero, not one. 2537 2538 @*/ 2539 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2540 { 2541 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2542 2543 PetscFunctionBegin; 2544 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2545 PetscValidPointer(indices,2); 2546 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2547 if (f) { 2548 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2549 } else { 2550 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2551 } 2552 PetscFunctionReturn(0); 2553 } 2554 2555 /* ----------------------------------------------------------------------------------------*/ 2556 2557 EXTERN_C_BEGIN 2558 #undef __FUNCT__ 2559 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2560 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2561 { 2562 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2563 PetscErrorCode ierr; 2564 size_t nz = aij->i[mat->rmap.n]; 2565 2566 PetscFunctionBegin; 2567 if (aij->nonew != 1) { 2568 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2569 } 2570 2571 /* allocate space for values if not already there */ 2572 if (!aij->saved_values) { 2573 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2574 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2575 } 2576 2577 /* copy values over */ 2578 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2579 PetscFunctionReturn(0); 2580 } 2581 EXTERN_C_END 2582 2583 #undef __FUNCT__ 2584 #define __FUNCT__ "MatStoreValues" 2585 /*@ 2586 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2587 example, reuse of the linear part of a Jacobian, while recomputing the 2588 nonlinear portion. 2589 2590 Collect on Mat 2591 2592 Input Parameters: 2593 . mat - the matrix (currently only AIJ matrices support this option) 2594 2595 Level: advanced 2596 2597 Common Usage, with SNESSolve(): 2598 $ Create Jacobian matrix 2599 $ Set linear terms into matrix 2600 $ Apply boundary conditions to matrix, at this time matrix must have 2601 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2602 $ boundary conditions again will not change the nonzero structure 2603 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2604 $ ierr = MatStoreValues(mat); 2605 $ Call SNESSetJacobian() with matrix 2606 $ In your Jacobian routine 2607 $ ierr = MatRetrieveValues(mat); 2608 $ Set nonlinear terms in matrix 2609 2610 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2611 $ // build linear portion of Jacobian 2612 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2613 $ ierr = MatStoreValues(mat); 2614 $ loop over nonlinear iterations 2615 $ ierr = MatRetrieveValues(mat); 2616 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2617 $ // call MatAssemblyBegin/End() on matrix 2618 $ Solve linear system with Jacobian 2619 $ endloop 2620 2621 Notes: 2622 Matrix must already be assemblied before calling this routine 2623 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2624 calling this routine. 2625 2626 When this is called multiple times it overwrites the previous set of stored values 2627 and does not allocated additional space. 2628 2629 .seealso: MatRetrieveValues() 2630 2631 @*/ 2632 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2633 { 2634 PetscErrorCode ierr,(*f)(Mat); 2635 2636 PetscFunctionBegin; 2637 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2638 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2639 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2640 2641 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2642 if (f) { 2643 ierr = (*f)(mat);CHKERRQ(ierr); 2644 } else { 2645 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2646 } 2647 PetscFunctionReturn(0); 2648 } 2649 2650 EXTERN_C_BEGIN 2651 #undef __FUNCT__ 2652 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2653 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2654 { 2655 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2656 PetscErrorCode ierr; 2657 PetscInt nz = aij->i[mat->rmap.n]; 2658 2659 PetscFunctionBegin; 2660 if (aij->nonew != 1) { 2661 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2662 } 2663 if (!aij->saved_values) { 2664 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2665 } 2666 /* copy values over */ 2667 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2668 PetscFunctionReturn(0); 2669 } 2670 EXTERN_C_END 2671 2672 #undef __FUNCT__ 2673 #define __FUNCT__ "MatRetrieveValues" 2674 /*@ 2675 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2676 example, reuse of the linear part of a Jacobian, while recomputing the 2677 nonlinear portion. 2678 2679 Collect on Mat 2680 2681 Input Parameters: 2682 . mat - the matrix (currently on AIJ matrices support this option) 2683 2684 Level: advanced 2685 2686 .seealso: MatStoreValues() 2687 2688 @*/ 2689 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2690 { 2691 PetscErrorCode ierr,(*f)(Mat); 2692 2693 PetscFunctionBegin; 2694 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2695 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2696 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2697 2698 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2699 if (f) { 2700 ierr = (*f)(mat);CHKERRQ(ierr); 2701 } else { 2702 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2703 } 2704 PetscFunctionReturn(0); 2705 } 2706 2707 2708 /* --------------------------------------------------------------------------------*/ 2709 #undef __FUNCT__ 2710 #define __FUNCT__ "MatCreateSeqAIJ" 2711 /*@C 2712 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2713 (the default parallel PETSc format). For good matrix assembly performance 2714 the user should preallocate the matrix storage by setting the parameter nz 2715 (or the array nnz). By setting these parameters accurately, performance 2716 during matrix assembly can be increased by more than a factor of 50. 2717 2718 Collective on MPI_Comm 2719 2720 Input Parameters: 2721 + comm - MPI communicator, set to PETSC_COMM_SELF 2722 . m - number of rows 2723 . n - number of columns 2724 . nz - number of nonzeros per row (same for all rows) 2725 - nnz - array containing the number of nonzeros in the various rows 2726 (possibly different for each row) or PETSC_NULL 2727 2728 Output Parameter: 2729 . A - the matrix 2730 2731 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2732 MatXXXXSetPreallocation() paradgm instead of this routine directly. This is definitely 2733 true if you plan to use the external direct solvers such as SuperLU, MUMPS or Spooles. 2734 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2735 2736 Notes: 2737 If nnz is given then nz is ignored 2738 2739 The AIJ format (also called the Yale sparse matrix format or 2740 compressed row storage), is fully compatible with standard Fortran 77 2741 storage. That is, the stored row and column indices can begin at 2742 either one (as in Fortran) or zero. See the users' manual for details. 2743 2744 Specify the preallocated storage with either nz or nnz (not both). 2745 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2746 allocation. For large problems you MUST preallocate memory or you 2747 will get TERRIBLE performance, see the users' manual chapter on matrices. 2748 2749 By default, this format uses inodes (identical nodes) when possible, to 2750 improve numerical efficiency of matrix-vector products and solves. We 2751 search for consecutive rows with the same nonzero structure, thereby 2752 reusing matrix information to achieve increased efficiency. 2753 2754 Options Database Keys: 2755 + -mat_no_inode - Do not use inodes 2756 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2757 - -mat_aij_oneindex - Internally use indexing starting at 1 2758 rather than 0. Note that when calling MatSetValues(), 2759 the user still MUST index entries starting at 0! 2760 2761 Level: intermediate 2762 2763 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2764 2765 @*/ 2766 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2767 { 2768 PetscErrorCode ierr; 2769 2770 PetscFunctionBegin; 2771 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2772 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2773 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2774 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2775 PetscFunctionReturn(0); 2776 } 2777 2778 #undef __FUNCT__ 2779 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2780 /*@C 2781 MatSeqAIJSetPreallocation - For good matrix assembly performance 2782 the user should preallocate the matrix storage by setting the parameter nz 2783 (or the array nnz). By setting these parameters accurately, performance 2784 during matrix assembly can be increased by more than a factor of 50. 2785 2786 Collective on MPI_Comm 2787 2788 Input Parameters: 2789 + B - The matrix 2790 . nz - number of nonzeros per row (same for all rows) 2791 - nnz - array containing the number of nonzeros in the various rows 2792 (possibly different for each row) or PETSC_NULL 2793 2794 Notes: 2795 If nnz is given then nz is ignored 2796 2797 The AIJ format (also called the Yale sparse matrix format or 2798 compressed row storage), is fully compatible with standard Fortran 77 2799 storage. That is, the stored row and column indices can begin at 2800 either one (as in Fortran) or zero. See the users' manual for details. 2801 2802 Specify the preallocated storage with either nz or nnz (not both). 2803 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2804 allocation. For large problems you MUST preallocate memory or you 2805 will get TERRIBLE performance, see the users' manual chapter on matrices. 2806 2807 You can call MatGetInfo() to get information on how effective the preallocation was; 2808 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2809 You can also run with the option -info and look for messages with the string 2810 malloc in them to see if additional memory allocation was needed. 2811 2812 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2813 entries or columns indices 2814 2815 By default, this format uses inodes (identical nodes) when possible, to 2816 improve numerical efficiency of matrix-vector products and solves. We 2817 search for consecutive rows with the same nonzero structure, thereby 2818 reusing matrix information to achieve increased efficiency. 2819 2820 Options Database Keys: 2821 + -mat_no_inode - Do not use inodes 2822 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2823 - -mat_aij_oneindex - Internally use indexing starting at 1 2824 rather than 0. Note that when calling MatSetValues(), 2825 the user still MUST index entries starting at 0! 2826 2827 Level: intermediate 2828 2829 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2830 2831 @*/ 2832 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2833 { 2834 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2835 2836 PetscFunctionBegin; 2837 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2838 if (f) { 2839 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2840 } 2841 PetscFunctionReturn(0); 2842 } 2843 2844 EXTERN_C_BEGIN 2845 #undef __FUNCT__ 2846 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2847 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2848 { 2849 Mat_SeqAIJ *b; 2850 PetscTruth skipallocation = PETSC_FALSE; 2851 PetscErrorCode ierr; 2852 PetscInt i; 2853 2854 PetscFunctionBegin; 2855 2856 if (nz == MAT_SKIP_ALLOCATION) { 2857 skipallocation = PETSC_TRUE; 2858 nz = 0; 2859 } 2860 2861 B->rmap.bs = B->cmap.bs = 1; 2862 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 2863 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 2864 2865 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2866 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2867 if (nnz) { 2868 for (i=0; i<B->rmap.n; i++) { 2869 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2870 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); 2871 } 2872 } 2873 2874 B->preallocated = PETSC_TRUE; 2875 b = (Mat_SeqAIJ*)B->data; 2876 2877 if (!skipallocation) { 2878 if (!b->imax) { 2879 ierr = PetscMalloc2(B->rmap.n,PetscInt,&b->imax,B->rmap.n,PetscInt,&b->ilen);CHKERRQ(ierr); 2880 ierr = PetscLogObjectMemory(B,2*B->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 2881 } 2882 if (!nnz) { 2883 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2884 else if (nz <= 0) nz = 1; 2885 for (i=0; i<B->rmap.n; i++) b->imax[i] = nz; 2886 nz = nz*B->rmap.n; 2887 } else { 2888 nz = 0; 2889 for (i=0; i<B->rmap.n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2890 } 2891 /* b->ilen will count nonzeros in each row so far. */ 2892 for (i=0; i<B->rmap.n; i++) { b->ilen[i] = 0; } 2893 2894 /* allocate the matrix space */ 2895 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 2896 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.n+1,PetscInt,&b->i);CHKERRQ(ierr); 2897 ierr = PetscLogObjectMemory(B,(B->rmap.n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 2898 b->i[0] = 0; 2899 for (i=1; i<B->rmap.n+1; i++) { 2900 b->i[i] = b->i[i-1] + b->imax[i-1]; 2901 } 2902 b->singlemalloc = PETSC_TRUE; 2903 b->free_a = PETSC_TRUE; 2904 b->free_ij = PETSC_TRUE; 2905 } else { 2906 b->free_a = PETSC_FALSE; 2907 b->free_ij = PETSC_FALSE; 2908 } 2909 2910 b->nz = 0; 2911 b->maxnz = nz; 2912 B->info.nz_unneeded = (double)b->maxnz; 2913 PetscFunctionReturn(0); 2914 } 2915 EXTERN_C_END 2916 2917 #undef __FUNCT__ 2918 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2919 /*@ 2920 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2921 2922 Input Parameters: 2923 + B - the matrix 2924 . i - the indices into j for the start of each row (starts with zero) 2925 . j - the column indices for each row (starts with zero) these must be sorted for each row 2926 - v - optional values in the matrix 2927 2928 Contributed by: Lisandro Dalchin 2929 2930 Level: developer 2931 2932 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 2933 2934 .keywords: matrix, aij, compressed row, sparse, sequential 2935 2936 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2937 @*/ 2938 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2939 { 2940 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2941 PetscErrorCode ierr; 2942 2943 PetscFunctionBegin; 2944 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2945 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2946 if (f) { 2947 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2948 } 2949 PetscFunctionReturn(0); 2950 } 2951 2952 EXTERN_C_BEGIN 2953 #undef __FUNCT__ 2954 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2955 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 2956 { 2957 PetscInt i; 2958 PetscInt m,n; 2959 PetscInt nz; 2960 PetscInt *nnz, nz_max = 0; 2961 PetscScalar *values; 2962 PetscErrorCode ierr; 2963 2964 PetscFunctionBegin; 2965 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2966 2967 if (Ii[0]) { 2968 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 2969 } 2970 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2971 for(i = 0; i < m; i++) { 2972 nz = Ii[i+1]- Ii[i]; 2973 nz_max = PetscMax(nz_max, nz); 2974 if (nz < 0) { 2975 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2976 } 2977 nnz[i] = nz; 2978 } 2979 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2980 ierr = PetscFree(nnz);CHKERRQ(ierr); 2981 2982 if (v) { 2983 values = (PetscScalar*) v; 2984 } else { 2985 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2986 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2987 } 2988 2989 for(i = 0; i < m; i++) { 2990 nz = Ii[i+1] - Ii[i]; 2991 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2992 } 2993 2994 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2995 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2996 2997 if (!v) { 2998 ierr = PetscFree(values);CHKERRQ(ierr); 2999 } 3000 PetscFunctionReturn(0); 3001 } 3002 EXTERN_C_END 3003 3004 /*MC 3005 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3006 based on compressed sparse row format. 3007 3008 Options Database Keys: 3009 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3010 3011 Level: beginner 3012 3013 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3014 M*/ 3015 3016 EXTERN_C_BEGIN 3017 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3018 EXTERN_C_END 3019 3020 EXTERN_C_BEGIN 3021 #undef __FUNCT__ 3022 #define __FUNCT__ "MatCreate_SeqAIJ" 3023 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3024 { 3025 Mat_SeqAIJ *b; 3026 PetscErrorCode ierr; 3027 PetscMPIInt size; 3028 3029 PetscFunctionBegin; 3030 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3031 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3032 3033 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3034 B->data = (void*)b; 3035 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3036 B->factor = 0; 3037 B->mapping = 0; 3038 b->row = 0; 3039 b->col = 0; 3040 b->icol = 0; 3041 b->reallocs = 0; 3042 b->ignorezeroentries = PETSC_FALSE; 3043 b->roworiented = PETSC_TRUE; 3044 b->nonew = 0; 3045 b->diag = 0; 3046 b->solve_work = 0; 3047 B->spptr = 0; 3048 b->saved_values = 0; 3049 b->idiag = 0; 3050 b->mdiag = 0; 3051 b->ssor_work = 0; 3052 b->omega = 1.0; 3053 b->fshift = 0.0; 3054 b->idiagvalid = PETSC_FALSE; 3055 b->keepzeroedrows = PETSC_FALSE; 3056 b->xtoy = 0; 3057 b->XtoY = 0; 3058 b->compressedrow.use = PETSC_FALSE; 3059 b->compressedrow.nrows = B->rmap.n; 3060 b->compressedrow.i = PETSC_NULL; 3061 b->compressedrow.rindex = PETSC_NULL; 3062 b->compressedrow.checked = PETSC_FALSE; 3063 B->same_nonzero = PETSC_FALSE; 3064 3065 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3066 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3067 "MatSeqAIJSetColumnIndices_SeqAIJ", 3068 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3069 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3070 "MatStoreValues_SeqAIJ", 3071 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3072 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3073 "MatRetrieveValues_SeqAIJ", 3074 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3075 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3076 "MatConvert_SeqAIJ_SeqSBAIJ", 3077 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3078 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3079 "MatConvert_SeqAIJ_SeqBAIJ", 3080 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3081 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3082 "MatConvert_SeqAIJ_SeqCSRPERM", 3083 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3084 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3085 "MatConvert_SeqAIJ_SeqCRL", 3086 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3087 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3088 "MatIsTranspose_SeqAIJ", 3089 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3090 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3091 "MatIsHermitianTranspose_SeqAIJ", 3092 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3093 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3094 "MatSeqAIJSetPreallocation_SeqAIJ", 3095 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3096 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3097 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3098 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3099 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3100 "MatReorderForNonzeroDiagonal_SeqAIJ", 3101 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3102 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 3103 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3104 PetscFunctionReturn(0); 3105 } 3106 EXTERN_C_END 3107 3108 #undef __FUNCT__ 3109 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3110 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3111 { 3112 Mat C; 3113 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3114 PetscErrorCode ierr; 3115 PetscInt i,m = A->rmap.n; 3116 3117 PetscFunctionBegin; 3118 *B = 0; 3119 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 3120 ierr = MatSetSizes(C,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 3121 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 3122 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3123 3124 ierr = PetscMapCopy(((PetscObject)A)->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 3125 ierr = PetscMapCopy(((PetscObject)A)->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 3126 3127 c = (Mat_SeqAIJ*)C->data; 3128 3129 C->factor = A->factor; 3130 3131 c->row = 0; 3132 c->col = 0; 3133 c->icol = 0; 3134 c->reallocs = 0; 3135 3136 C->assembled = PETSC_TRUE; 3137 3138 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3139 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3140 for (i=0; i<m; i++) { 3141 c->imax[i] = a->imax[i]; 3142 c->ilen[i] = a->ilen[i]; 3143 } 3144 3145 /* allocate the matrix space */ 3146 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3147 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3148 c->singlemalloc = PETSC_TRUE; 3149 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3150 if (m > 0) { 3151 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3152 if (cpvalues == MAT_COPY_VALUES) { 3153 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3154 } else { 3155 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3156 } 3157 } 3158 3159 c->ignorezeroentries = a->ignorezeroentries; 3160 c->roworiented = a->roworiented; 3161 c->nonew = a->nonew; 3162 if (a->diag) { 3163 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3164 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3165 for (i=0; i<m; i++) { 3166 c->diag[i] = a->diag[i]; 3167 } 3168 } else c->diag = 0; 3169 c->solve_work = 0; 3170 c->saved_values = 0; 3171 c->idiag = 0; 3172 c->ssor_work = 0; 3173 c->keepzeroedrows = a->keepzeroedrows; 3174 c->free_a = PETSC_TRUE; 3175 c->free_ij = PETSC_TRUE; 3176 c->xtoy = 0; 3177 c->XtoY = 0; 3178 3179 c->nz = a->nz; 3180 c->maxnz = a->maxnz; 3181 C->preallocated = PETSC_TRUE; 3182 3183 c->compressedrow.use = a->compressedrow.use; 3184 c->compressedrow.nrows = a->compressedrow.nrows; 3185 c->compressedrow.checked = a->compressedrow.checked; 3186 if ( a->compressedrow.checked && a->compressedrow.use){ 3187 i = a->compressedrow.nrows; 3188 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3189 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3190 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3191 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3192 } else { 3193 c->compressedrow.use = PETSC_FALSE; 3194 c->compressedrow.i = PETSC_NULL; 3195 c->compressedrow.rindex = PETSC_NULL; 3196 } 3197 C->same_nonzero = A->same_nonzero; 3198 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3199 3200 *B = C; 3201 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3202 PetscFunctionReturn(0); 3203 } 3204 3205 #undef __FUNCT__ 3206 #define __FUNCT__ "MatLoad_SeqAIJ" 3207 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 3208 { 3209 Mat_SeqAIJ *a; 3210 Mat B; 3211 PetscErrorCode ierr; 3212 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3213 int fd; 3214 PetscMPIInt size; 3215 MPI_Comm comm; 3216 3217 PetscFunctionBegin; 3218 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3219 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3220 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3221 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3222 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3223 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3224 M = header[1]; N = header[2]; nz = header[3]; 3225 3226 if (nz < 0) { 3227 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3228 } 3229 3230 /* read in row lengths */ 3231 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3232 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3233 3234 /* check if sum of rowlengths is same as nz */ 3235 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3236 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3237 3238 /* create our matrix */ 3239 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3240 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3241 ierr = MatSetType(B,type);CHKERRQ(ierr); 3242 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3243 a = (Mat_SeqAIJ*)B->data; 3244 3245 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3246 3247 /* read in nonzero values */ 3248 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3249 3250 /* set matrix "i" values */ 3251 a->i[0] = 0; 3252 for (i=1; i<= M; i++) { 3253 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3254 a->ilen[i-1] = rowlengths[i-1]; 3255 } 3256 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3257 3258 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3259 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3260 *A = B; 3261 PetscFunctionReturn(0); 3262 } 3263 3264 #undef __FUNCT__ 3265 #define __FUNCT__ "MatEqual_SeqAIJ" 3266 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3267 { 3268 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3269 PetscErrorCode ierr; 3270 3271 PetscFunctionBegin; 3272 /* If the matrix dimensions are not equal,or no of nonzeros */ 3273 if ((A->rmap.n != B->rmap.n) || (A->cmap.n != B->cmap.n) ||(a->nz != b->nz)) { 3274 *flg = PETSC_FALSE; 3275 PetscFunctionReturn(0); 3276 } 3277 3278 /* if the a->i are the same */ 3279 ierr = PetscMemcmp(a->i,b->i,(A->rmap.n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3280 if (!*flg) PetscFunctionReturn(0); 3281 3282 /* if a->j are the same */ 3283 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3284 if (!*flg) PetscFunctionReturn(0); 3285 3286 /* if a->a are the same */ 3287 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3288 3289 PetscFunctionReturn(0); 3290 3291 } 3292 3293 #undef __FUNCT__ 3294 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3295 /*@ 3296 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3297 provided by the user. 3298 3299 Collective on MPI_Comm 3300 3301 Input Parameters: 3302 + comm - must be an MPI communicator of size 1 3303 . m - number of rows 3304 . n - number of columns 3305 . i - row indices 3306 . j - column indices 3307 - a - matrix values 3308 3309 Output Parameter: 3310 . mat - the matrix 3311 3312 Level: intermediate 3313 3314 Notes: 3315 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3316 once the matrix is destroyed 3317 3318 You cannot set new nonzero locations into this matrix, that will generate an error. 3319 3320 The i and j indices are 0 based 3321 3322 The format which is used for the sparse matrix input, is equivalent to a 3323 row-major ordering.. i.e for the following matrix, the input data expected is 3324 as shown: 3325 3326 1 0 0 3327 2 0 3 3328 4 5 6 3329 3330 i = {0,1,3,6} [size = nrow+1 = 3+1] 3331 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3332 v = {1,2,3,4,5,6} [size = nz = 6] 3333 3334 3335 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3336 3337 @*/ 3338 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3339 { 3340 PetscErrorCode ierr; 3341 PetscInt ii; 3342 Mat_SeqAIJ *aij; 3343 #if defined(PETSC_USE_DEBUG) 3344 PetscInt jj; 3345 #endif 3346 3347 PetscFunctionBegin; 3348 if (i[0]) { 3349 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3350 } 3351 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3352 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3353 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3354 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3355 aij = (Mat_SeqAIJ*)(*mat)->data; 3356 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3357 3358 aij->i = i; 3359 aij->j = j; 3360 aij->a = a; 3361 aij->singlemalloc = PETSC_FALSE; 3362 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3363 aij->free_a = PETSC_FALSE; 3364 aij->free_ij = PETSC_FALSE; 3365 3366 for (ii=0; ii<m; ii++) { 3367 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3368 #if defined(PETSC_USE_DEBUG) 3369 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]); 3370 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3371 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); 3372 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); 3373 } 3374 #endif 3375 } 3376 #if defined(PETSC_USE_DEBUG) 3377 for (ii=0; ii<aij->i[m]; ii++) { 3378 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3379 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3380 } 3381 #endif 3382 3383 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3384 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3385 PetscFunctionReturn(0); 3386 } 3387 3388 #undef __FUNCT__ 3389 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3390 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3391 { 3392 PetscErrorCode ierr; 3393 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3394 3395 PetscFunctionBegin; 3396 if (coloring->ctype == IS_COLORING_GLOBAL) { 3397 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3398 a->coloring = coloring; 3399 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3400 PetscInt i,*larray; 3401 ISColoring ocoloring; 3402 ISColoringValue *colors; 3403 3404 /* set coloring for diagonal portion */ 3405 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3406 for (i=0; i<A->cmap.n; i++) { 3407 larray[i] = i; 3408 } 3409 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap.n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3410 ierr = PetscMalloc((A->cmap.n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3411 for (i=0; i<A->cmap.n; i++) { 3412 colors[i] = coloring->colors[larray[i]]; 3413 } 3414 ierr = PetscFree(larray);CHKERRQ(ierr); 3415 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap.n,colors,&ocoloring);CHKERRQ(ierr); 3416 a->coloring = ocoloring; 3417 } 3418 PetscFunctionReturn(0); 3419 } 3420 3421 #if defined(PETSC_HAVE_ADIC) 3422 EXTERN_C_BEGIN 3423 #include "adic/ad_utils.h" 3424 EXTERN_C_END 3425 3426 #undef __FUNCT__ 3427 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3428 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3429 { 3430 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3431 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3432 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3433 ISColoringValue *color; 3434 3435 PetscFunctionBegin; 3436 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3437 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3438 color = a->coloring->colors; 3439 /* loop over rows */ 3440 for (i=0; i<m; i++) { 3441 nz = ii[i+1] - ii[i]; 3442 /* loop over columns putting computed value into matrix */ 3443 for (j=0; j<nz; j++) { 3444 *v++ = values[color[*jj++]]; 3445 } 3446 values += nlen; /* jump to next row of derivatives */ 3447 } 3448 PetscFunctionReturn(0); 3449 } 3450 #endif 3451 3452 #undef __FUNCT__ 3453 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3454 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3455 { 3456 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3457 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j; 3458 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3459 ISColoringValue *color; 3460 3461 PetscFunctionBegin; 3462 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3463 color = a->coloring->colors; 3464 /* loop over rows */ 3465 for (i=0; i<m; i++) { 3466 nz = ii[i+1] - ii[i]; 3467 /* loop over columns putting computed value into matrix */ 3468 for (j=0; j<nz; j++) { 3469 *v++ = values[color[*jj++]]; 3470 } 3471 values += nl; /* jump to next row of derivatives */ 3472 } 3473 PetscFunctionReturn(0); 3474 } 3475 3476 /* 3477 Special version for direct calls from Fortran 3478 */ 3479 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3480 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3481 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3482 #define matsetvaluesseqaij_ matsetvaluesseqaij 3483 #endif 3484 3485 /* Change these macros so can be used in void function */ 3486 #undef CHKERRQ 3487 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3488 #undef SETERRQ2 3489 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(((PetscObject)A)->comm,ierr) 3490 3491 EXTERN_C_BEGIN 3492 #undef __FUNCT__ 3493 #define __FUNCT__ "matsetvaluesseqaij_" 3494 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3495 { 3496 Mat A = *AA; 3497 PetscInt m = *mm, n = *nn; 3498 InsertMode is = *isis; 3499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3500 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3501 PetscInt *imax,*ai,*ailen; 3502 PetscErrorCode ierr; 3503 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3504 PetscScalar *ap,value,*aa; 3505 PetscTruth ignorezeroentries = a->ignorezeroentries; 3506 PetscTruth roworiented = a->roworiented; 3507 3508 PetscFunctionBegin; 3509 MatPreallocated(A); 3510 imax = a->imax; 3511 ai = a->i; 3512 ailen = a->ilen; 3513 aj = a->j; 3514 aa = a->a; 3515 3516 for (k=0; k<m; k++) { /* loop over added rows */ 3517 row = im[k]; 3518 if (row < 0) continue; 3519 #if defined(PETSC_USE_DEBUG) 3520 if (row >= A->rmap.n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3521 #endif 3522 rp = aj + ai[row]; ap = aa + ai[row]; 3523 rmax = imax[row]; nrow = ailen[row]; 3524 low = 0; 3525 high = nrow; 3526 for (l=0; l<n; l++) { /* loop over added columns */ 3527 if (in[l] < 0) continue; 3528 #if defined(PETSC_USE_DEBUG) 3529 if (in[l] >= A->cmap.n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3530 #endif 3531 col = in[l]; 3532 if (roworiented) { 3533 value = v[l + k*n]; 3534 } else { 3535 value = v[k + l*m]; 3536 } 3537 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3538 3539 if (col <= lastcol) low = 0; else high = nrow; 3540 lastcol = col; 3541 while (high-low > 5) { 3542 t = (low+high)/2; 3543 if (rp[t] > col) high = t; 3544 else low = t; 3545 } 3546 for (i=low; i<high; i++) { 3547 if (rp[i] > col) break; 3548 if (rp[i] == col) { 3549 if (is == ADD_VALUES) ap[i] += value; 3550 else ap[i] = value; 3551 goto noinsert; 3552 } 3553 } 3554 if (value == 0.0 && ignorezeroentries) goto noinsert; 3555 if (nonew == 1) goto noinsert; 3556 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3557 MatSeqXAIJReallocateAIJ(A,A->rmap.n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3558 N = nrow++ - 1; a->nz++; high++; 3559 /* shift up all the later entries in this row */ 3560 for (ii=N; ii>=i; ii--) { 3561 rp[ii+1] = rp[ii]; 3562 ap[ii+1] = ap[ii]; 3563 } 3564 rp[i] = col; 3565 ap[i] = value; 3566 noinsert:; 3567 low = i + 1; 3568 } 3569 ailen[row] = nrow; 3570 } 3571 A->same_nonzero = PETSC_FALSE; 3572 PetscFunctionReturnVoid(); 3573 } 3574 EXTERN_C_END 3575