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