- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi, I have used The intel MKL to calculate an eigenvalue equation for Self Consistent Field Approximation (Hartree Fock method). The program part is below. But For a degenerate matrix, Sometimes it changes order of eigenvectors like that
tempCC1: [ 0.999927 -0.012047 0.000000 -0.001017 0.000000 ]
tempCC2: [ 0.010611 0.914781 0.000000 -0.403812 0.000000 ]
tempCC3: [ 0.000000 0.000000 0.000000 0.000000 1.000000 ]
tempCC4: [ 0.005795 0.403771 0.000000 0.914842 0.000000 ]
tempCC5: [ 0.000000 0.000000 1.000000 0.000000 0.000000 ]
whereas it should make a progress like that
Eigenvalues:
en[1] = -11.301550274377
en[2] = -0.677494616682
en[3] = -0.003523316893
en[4] = -0.003523316893
en[5] = -0.003523316893
Eigenvectors (rows):
tempCC1: [ 0.999938 -0.011161 0.000000 0.000000 0.000000 ]
tempCC2: [ 0.011161 0.999938 0.000000 0.000000 0.000000 ]
tempCC3: [ 0.000000 0.000000 1.000000 0.000000 0.000000 ]
tempCC4: [ 0.000000 0.000000 0.000000 1.000000 0.000000 ]
tempCC5: [ 0.000000 0.000000 0.000000 0.000000 1.000000 ]
I want to know whether or not I make a mistake in the program part below. I hope someone helps me and corrects it. thank you in advance, SEDSIL.
const
LAPACK_COL_MAJOR = 102;
MKL_PATH = 'pathto\Intel\oneAPI\mkl\latest\bin\mkl_rt.dll';
OMP_PATH = 'pathto\Intel\oneAPI\compiler\latest\bin\libiomp5md.dll';
type
TLAPACKE_dsyev = function(layout: Integer; jobz: Char; uplo: Char;
n: Integer; a: PDouble; lda: Integer;
w: PDouble): Integer; cdecl;
function CalculateEigens(const a: RealArrayNPbyNP;
var d: RealArrayNP;
var v: RealArrayNPbyNP): Boolean;
var
i, j, ret: Integer;
a_data: array of Double;
w: array of Double;
hMKL, hOMP: THandle;
f_dsyev: TLAPACKE_dsyev;
begin
Result := False;
hMKL := 0;
hOMP := 0;
try
SetLength(a_data, np*np);
SetLength(w, np);
for i := 1 to np do
for j := 1 to np do
a_data[(j-1)*np + (i-1)] := a[i,j];
hOMP := LoadLibrary(PChar(OMP_PATH));
hMKL := LoadLibrary(PChar(MKL_PATH));
if (hOMP = 0) or (hMKL = 0) then Exit;
f_dsyev := TLAPACKE_dsyev(GetProcAddress(hMKL, 'LAPACKE_dsyev'));
if not Assigned(f_dsyev) then Exit;
ret := f_dsyev(LAPACK_COL_MAJOR, 'V', 'U', np, @a_data[0], np, @w[0]);
if ret <> 0 then Exit;
for i := 1 to np do
begin
d[i] := w[i-1];
for j := 1 to np do
v[i,j] := a_data[(j-1)*np + (i-1)];
end;
Result := True;
finally
if hMKL <> 0 then FreeLibrary(hMKL);
if hOMP <> 0 then FreeLibrary(hOMP);
end;
end;
Link Copied
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page