| EXTERN PetscErrorCode IPInnerProduct(IP ip,Vec,Vec,PetscScalar*); |
| EXTERN PetscErrorCode IPInnerProductBegin(IP ip,Vec,Vec,PetscScalar*); |
| EXTERN PetscErrorCode IPInnerProductEnd(IP ip,Vec,Vec,PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProduct(IP ip,PetscInt,Vec,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProductBegin(IP ip,PetscInt,Vec,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProductEnd(IP ip,PetscInt,Vec,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProduct(IP ip,Vec,PetscInt,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProductBegin(IP ip,Vec,PetscInt,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPMInnerProductEnd(IP ip,Vec,PetscInt,const Vec[],PetscScalar*); |
| EXTERN PetscErrorCode IPNorm(IP ip,Vec,PetscReal*); |
| EXTERN PetscErrorCode IPNormBegin(IP ip,Vec,PetscReal*); |
| EXTERN PetscErrorCode IPNormEnd(IP ip,Vec,PetscReal*); |
| ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr); |
| ierr = IPOrthogonalize(eps->ip,eps->nds,PETSC_NULL,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL,eps->work[0]);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(eps->ip,j+1,f,V,H+m*j);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(eps->ip,f,j+1,V,H+m*j);CHKERRQ(ierr); |
| if (j>k) { |
| ierr = IPMInnerProductBegin(eps->ip,j,V[j],V,lhh);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(eps->ip,V[j],j,V,lhh);CHKERRQ(ierr); |
| ierr = IPInnerProductBegin(eps->ip,V[j],V[j],&dot);CHKERRQ(ierr); |
| } |
| if (j>k+1) { |
| ierr = VecDotBegin(u,V[j-2],&dot2);CHKERRQ(ierr); |
| } |
| ierr = IPMInnerProductEnd(eps->ip,j+1,f,V,H+m*j);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(eps->ip,f,j+1,V,H+m*j);CHKERRQ(ierr); |
| if (j>k) { |
| ierr = IPMInnerProductEnd(eps->ip,j,V[j],V,lhh);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(eps->ip,V[j],j,V,lhh);CHKERRQ(ierr); |
| ierr = IPInnerProductEnd(eps->ip,V[j],V[j],&dot);CHKERRQ(ierr); |
| } |
| if (j>k+1) { |
| ierr = VecCopy(t,V[m-1]);CHKERRQ(ierr); |
| H[m*(m-2)+m-1] = norm2; |
| ierr = IPMInnerProduct(eps->ip,m,f,V,lhh);CHKERRQ(ierr); |
| ierr = IPMInnerProduct(eps->ip,f,m,V,lhh);CHKERRQ(ierr); |
| ierr = VecSet(w,0.0);CHKERRQ(ierr); |
| ierr = VecMAXPY(w,m,lhh,V);CHKERRQ(ierr); |
| ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr); |
| ierr = IPOrthogonalize(eps->ip,eps->nds,PETSC_NULL,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL,eps->work[0]);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(eps->ip,j+1,f,V,H+m*j);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(eps->ip,f,j+1,V,H+m*j);CHKERRQ(ierr); |
| if (j>k) { |
| ierr = IPInnerProductBegin(eps->ip,V[j],V[j],&dot);CHKERRQ(ierr); |
| } |
| ierr = IPMInnerProductEnd(eps->ip,j+1,f,V,H+m*j);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(eps->ip,f,j+1,V,H+m*j);CHKERRQ(ierr); |
| if (j>k) { |
| ierr = IPInnerProductEnd(eps->ip,V[j],V[j],&dot);CHKERRQ(ierr); |
| } |
| Input Parameters: |
| + ip - the spectral transformation context |
| . x - the first input vector |
| . n - number of vectors in y |
| . x - the first input vector |
| - y - array of vectors |
| Output Parameter: |
| .seealso: IPSetBilinearForm(), IPApplyB(), VecMDot(), IPInnerProduct() |
| @*/ |
| PetscErrorCode IPMInnerProduct(IP ip,PetscInt n,Vec x,const Vec y[],PetscScalar *p) |
| PetscErrorCode IPMInnerProduct(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p) |
| { |
| PetscErrorCode ierr; |
| Input Parameters: |
| + ip - the spectral transformation context |
| . x - the first input vector |
| . n - number of vectors in y |
| . x - the first input vector |
| . y - array of vectors |
| - p - where the result will go |
| IPNormEnd(), IPInnerProduct() |
| @*/ |
| PetscErrorCode IPMInnerProductBegin(IP ip,PetscInt n,Vec x,const Vec y[],PetscScalar *p) |
| PetscErrorCode IPMInnerProductBegin(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p) |
| { |
| PetscErrorCode ierr; |
| Input Parameters: |
| + ip - the spectral transformation context |
| . x - the first input vector |
| . n - number of vectors in y |
| . x - the first input vector |
| - y - array of vectors |
| Output Parameter: |
| IPNormEnd(), IPInnerProduct() |
| @*/ |
| PetscErrorCode IPMInnerProductEnd(IP ip,PetscInt n,Vec x,const Vec y[],PetscScalar *p) |
| PetscErrorCode IPMInnerProductEnd(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p) |
| { |
| PetscErrorCode ierr; |
| if (onorm || norm) { ierr = IPInnerProductEnd(ip,v,v,&alpha);CHKERRQ(ierr); } |
| } else { /* merge comunications */ |
| if (onorm || norm) { |
| ierr = IPMInnerProductBegin(ip,n,v,V,H);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(ip,v,n,V,H);CHKERRQ(ierr); |
| ierr = IPInnerProductBegin(ip,v,v,&alpha);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(ip,n,v,V,H);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(ip,v,n,V,H);CHKERRQ(ierr); |
| ierr = IPInnerProductEnd(ip,v,v,&alpha);CHKERRQ(ierr); |
| } else { /* use simpler function */ |
| ierr = IPMInnerProduct(ip,n,v,V,H);CHKERRQ(ierr); |
| ierr = IPMInnerProduct(ip,v,n,V,H);CHKERRQ(ierr); |
| } |
| } |
| + ip - the inner product (IP) context |
| . n - number of columns of V |
| . which - logical array indicating columns of V to be used |
| - V - set of vectors |
| . V - set of vectors |
| - work - workspace vector |
| Input/Output Parameter: |
| . v - (input) vector to be orthogonalized and (output) result of |
| . V - set of vectors |
| . m - starting column |
| . n - ending column |
| - ldr - leading dimension of R |
| . ldr - leading dimension of R |
| - work - workspace vector |
| Output Parameter: |
| . R - triangular matrix of QR factorization |
| ierr = VecDuplicate(v,&w);CHKERRQ(ierr); |
| for (j=0;j<n;j++) { |
| ierr = IPMInnerProduct(ip,n,V[j],W,vw+j*n);CHKERRQ(ierr); |
| ierr = IPMInnerProduct(ip,V[j],n,W,vw+j*n);CHKERRQ(ierr); |
| } |
| lwork = n; |
| ierr = PetscMalloc(n*sizeof(PetscScalar),&tau);CHKERRQ(ierr); |
| /* h = W^* v */ |
| /* q = v - V h */ |
| ierr = IPMInnerProduct(ip,n,v,W,H);CHKERRQ(ierr); |
| ierr = IPMInnerProduct(ip,v,n,W,H);CHKERRQ(ierr); |
| BLAStrsm_("L","L","N","N",&n,&ione,&one,vw,&n,H,&n,1,1,1,1); |
| LAPACKormlq_("L","N",&n,&ione,&n,vw,&n,tau,H,&n,work,&lwork,&info,1,1); |
| if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORMLQ %i",info); |
| for (i=k+1;i<n;i++) { |
| ierr = SVDMatMult(svd,PETSC_TRUE,U[i-1],V[i]);CHKERRQ(ierr); |
| ierr = IPNormBegin(svd->ip,U[i-1],&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,i,V[i],V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,V[i],i,V,work);CHKERRQ(ierr); |
| ierr = IPNormEnd(svd->ip,U[i-1],&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,i,V[i],V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,V[i],i,V,work);CHKERRQ(ierr); |
| ierr = VecScale(U[i-1],1.0/a);CHKERRQ(ierr); |
| ierr = VecScale(V[i],1.0/a);CHKERRQ(ierr); |
| } |
| ierr = SVDMatMult(svd,PETSC_TRUE,U[n-1],v);CHKERRQ(ierr); |
| ierr = IPNormBegin(svd->ip,U[n-1],&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,n,v,V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,v,n,V,work);CHKERRQ(ierr); |
| ierr = IPNormEnd(svd->ip,U[n-1],&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,n,v,V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,v,n,V,work);CHKERRQ(ierr); |
| ierr = VecScale(U[n-1],1.0/a);CHKERRQ(ierr); |
| ierr = VecScale(v,1.0/a);CHKERRQ(ierr); |
| for (i=k+1;i<n;i++) { |
| ierr = SVDMatMult(svd,PETSC_TRUE,u,V[i]);CHKERRQ(ierr); |
| ierr = IPNormBegin(svd->ip,u,&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,i,V[i],V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,V[i],i,V,work);CHKERRQ(ierr); |
| ierr = IPNormEnd(svd->ip,u,&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,i,V[i],V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,V[i],i,V,work);CHKERRQ(ierr); |
| ierr = VecScale(u,1.0/a);CHKERRQ(ierr); |
| ierr = VecScale(V[i],1.0/a);CHKERRQ(ierr); |
| } |
| ierr = SVDMatMult(svd,PETSC_TRUE,u,v);CHKERRQ(ierr); |
| ierr = IPNormBegin(svd->ip,u,&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,n,v,V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductBegin(svd->ip,v,n,V,work);CHKERRQ(ierr); |
| ierr = IPNormEnd(svd->ip,u,&a);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,n,v,V,work);CHKERRQ(ierr); |
| ierr = IPMInnerProductEnd(svd->ip,v,n,V,work);CHKERRQ(ierr); |
| ierr = VecScale(u,1.0/a);CHKERRQ(ierr); |
| ierr = VecScale(v,1.0/a);CHKERRQ(ierr); |