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