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