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