| Line 225... |
Line 225... |
LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info);
|
LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGETRI %d",info);
|
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGETRI %d",info);
|
hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
|
hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
|
|
|
*cond = hn*hin;
|
*cond = hn*hin;
|
|
PetscFunctionReturn(0);
|
|
#endif
|
|
}
|
|
|
|
#undef __FUNCT__
|
|
#define __FUNCT__ "PSTranslateHarmonic_NHEP"
|
|
PetscErrorCode PSTranslateHarmonic_NHEP(PS ps,PetscScalar tau,PetscScalar beta,PetscScalar *g)
|
|
{
|
|
#if defined(PETSC_MISSING_LAPACK_GETRF) || defined(PETSC_MISSING_LAPACK_GETRS)
|
|
PetscFunctionBegin;
|
|
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF/GETRS - Lapack routines are unavailable.");
|
|
#else
|
|
PetscErrorCode ierr;
|
|
PetscInt i;
|
|
PetscBLASInt *ipiv,info,n,ld,one=1;
|
|
PetscScalar *A,*B;
|
|
|
|
PetscFunctionBegin;
|
|
n = PetscBLASIntCast(ps->n);
|
|
ld = PetscBLASIntCast(ps->ld);
|
|
ierr = PSAllocateWork_Private(ps,0,0,ld);CHKERRQ(ierr);
|
|
ipiv = ps->iwork;
|
|
/* use workspace matrix W to factor A-tau*eye(n) */
|
|
if (!ps->mat[PS_MAT_W]) {
|
|
ierr = PSAllocateMat_Private(ps,PS_MAT_W);CHKERRQ(ierr);
|
|
}
|
|
B = ps->mat[PS_MAT_W];
|
|
A = ps->mat[PS_MAT_A];
|
|
ierr = PetscMemcpy(B,A,sizeof(PetscScalar)*ld*ld);CHKERRQ(ierr);
|
|
|
|
/* Vector g initialy stores b = beta*e_n^T */
|
|
ierr = PetscMemzero(g,n*sizeof(PetscScalar));CHKERRQ(ierr);
|
|
g[n-1] = beta;
|
|
|
|
/* g = (A-sigma*eye(n))'\b */
|
|
for (i=0;i<n;i++)
|
|
B[i+i*ld] -= tau;
|
|
LAPACKgetrf_(&n,&n,B,&ld,ipiv,&info);
|
|
if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
|
|
if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
|
|
ierr = PetscLogFlops(2.0*n*n*n/3.0);CHKERRQ(ierr);
|
|
LAPACKgetrs_("C",&n,&one,B,&ld,ipiv,g,&ld,&info);
|
|
if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
|
|
ierr = PetscLogFlops(2.0*n*n-n);CHKERRQ(ierr);
|
|
|
|
/* A = A + g*b' */
|
|
for (i=0;i<n;i++)
|
|
A[i+(n-1)*ld] = A[i+(n-1)*ld] + g[i]*beta;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
EXTERN_C_BEGIN
|
EXTERN_C_BEGIN
|
| Line 240... |
Line 288... |
ps->ops->allocate = PSAllocate_NHEP;
|
ps->ops->allocate = PSAllocate_NHEP;
|
//ps->ops->computevector = PSComputeVector_NHEP;
|
//ps->ops->computevector = PSComputeVector_NHEP;
|
ps->ops->solve = PSSolve_NHEP;
|
ps->ops->solve = PSSolve_NHEP;
|
ps->ops->sort = PSSort_NHEP;
|
ps->ops->sort = PSSort_NHEP;
|
ps->ops->cond = PSCond_NHEP;
|
ps->ops->cond = PSCond_NHEP;
|
|
ps->ops->translate = PSTranslateHarmonic_NHEP;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
EXTERN_C_END
|
EXTERN_C_END
|
|
|
|
|