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