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