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