Commit c688e23a by Francois Gygi

use LAPACK function dsyev

git-svn-id: http://qboxcode.org/svn/qb/trunk@1458 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 0fc8d558
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// D3vector.h
//
// double 3-vectors
//
////////////////////////////////////////////////////////////////////////////////
// $Id: D3vector.h,v 1.8 2008/09/08 15:56:18 fgygi Exp $
#ifndef D3VECTOR_H
#define D3VECTOR_H
#include <iostream>
#include <cmath>
#include <cassert>
class D3vector
{
public:
double x, y, z;
// explicit constructor to avoid implicit conversion from double to D3vector
explicit D3vector(const double& xv, const double& yv, const double& zv) :
x(xv), y(yv), z(zv) {}
explicit D3vector(void) : x(0.0), y(0.0), z(0.0) {}
explicit D3vector(const double* r) : x(r[0]), y(r[1]), z(r[2]) {}
double& operator[](const int &i)
{
assert(i>=0 && i <3);
if ( i == 0 ) return x;
else if ( i == 1 ) return y;
else return z;
}
bool operator==(const D3vector &rhs) const
{
return x == rhs.x && y == rhs.y && z == rhs.z;
}
bool operator!=(const D3vector &rhs) const
{
return x != rhs.x || y != rhs.y || z != rhs.z;
}
D3vector& operator += ( const D3vector& rhs )
{
x += rhs.x; y += rhs.y; z += rhs.z;
return *this;
}
D3vector& operator -= ( const D3vector& rhs )
{
x -= rhs.x; y -= rhs.y; z -= rhs.z;
return *this;
}
D3vector& operator *= ( const double& rhs )
{
x *= rhs; y *= rhs; z *= rhs;
return *this;
}
D3vector& operator /= ( const double& rhs )
{
x /= rhs; y /= rhs; z /= rhs;
return *this;
}
friend const D3vector operator + (const D3vector& lhs, const D3vector& rhs )
{
return D3vector(lhs) += rhs;
}
friend const D3vector operator - ( const D3vector& a, const D3vector& b )
{
return D3vector(a) -= b;
}
friend D3vector operator - ( const D3vector& a ) // unary minus
{
return D3vector( -a.x, -a.y, -a.z );
}
friend D3vector operator * ( const double& a, const D3vector& b )
{
return D3vector(b) *= a;
}
friend D3vector operator * ( const D3vector& a, const double& b )
{
return D3vector(a) *= b;
}
friend D3vector operator / ( const D3vector& a, const double& b )
{
return D3vector(a) /= b;
}
// scalar product
friend double operator * ( const D3vector& a, const D3vector& b )
{
return a.x * b.x + a.y * b.y + a.z * b.z ;
}
friend D3vector operator ^ ( const D3vector& a, const D3vector& b )
{
return D3vector( a.y * b.z - a.z * b.y ,
a.z * b.x - a.x * b.z ,
a.x * b.y - a.y * b.x );
}
friend D3vector rotate ( const D3vector& x, const D3vector& w )
{
if ( length(x) == 0.0 ) return x; // x has zero length
double theta = length( w ); // rotate by zero
if ( theta == 0.0 ) return x;
D3vector ew = normalized ( w );
D3vector v = w ^ x;
if ( length( v ) == 0.0 ) return x; // x is parallel to the rotation axis
v = normalized( v );
D3vector u = v ^ ew;
double p = x * u;
return (x*ew)*ew + p*cos(theta)*u + p*sin(theta)*v ;
}
friend double length( const D3vector& a )
{
return sqrt( a.x * a.x + a.y * a.y + a.z * a.z );
}
friend double norm( const D3vector& a )
{
return a.x * a.x + a.y * a.y + a.z * a.z;
}
friend D3vector normalized( const D3vector a )
{
return a / length( a );
}
friend std::ostream& operator << ( std::ostream& os, const D3vector& v )
{
os << v.x << " " << v.y << " " << v.z;
return os;
}
friend std::istream& operator >> ( std::istream& is, D3vector& v )
{
is >> v.x >> v.y >> v.z ;
return is;
}
};
#endif
#-------------------------------------------------------------------------------
#
# Copyright (c) 2008 The Regents of the University of California
# Copyright (c) 2008-2014 The Regents of the University of California
#
# This file is part of Qbox
#
......@@ -11,37 +11,15 @@
# or <http://www.gnu.org/licenses/>.
#
#-------------------------------------------------------------------------------
# $Id: Makefile,v 1.65 2009/09/08 05:40:50 fgygi Exp $
#------------------------------------------------------------------------------
#
include $(TARGET).mk
#------------------------------------------------------------------------------
SRCDIR = ../../src/
EXEC = dynmat
# use source files in SRCDIR but not objects, which may use MPI
vpath %.h $(SRCDIR)
vpath %.C $(SRCDIR)
OBJECTS = dynmat.o Matrix.o Context.o
CXXFLAGS += -DTARGET='"$(TARGET)"' -I$(SRCDIR)
$(EXEC): $(OBJECTS)
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
#------------------------------------------------------------------------------
# generate dependencies in makefile: use -Y to avoid library header files
# that are likely to be different on other platforms.
depend :
makedepend -Y -D$(PLT) *.[cCh]
#makedepend -D$(PLT) -I/usr/include/CC *.[cC]
EXEC=dynmat
LDFLAGS=-llapack
CXXFLAGS=-g
$(EXEC): dynmat.C
$(CXX) -o $@ $^ $(LDFLAGS)
#------------------------------------------------------------------------------
# Cleanup object files
clean :
rm -f *.o
rm -f $(EXEC) $(TESTEXECS)
#------------------------------------------------------------------------------
ctags :
etags -o tags *.[Ch]
rm -f $(EXEC)
#------------------------------------------------------------------------------
# DO NOT DELETE
Context.o: Context.h
dynmat.o: Context.h Matrix.h
Matrix.o: Context.h Matrix.h
......@@ -40,8 +40,9 @@
#include <vector>
using namespace std;
#include "Context.h"
#include "Matrix.h"
extern "C"
void dsyev_(const char *jobz, const char *uplo, const int *n, double *a,
const int *lda, double *w, double *wrk, int *lwrk, int *info);
int main(int argc, char **argv)
{
......@@ -95,9 +96,8 @@ int main(int argc, char **argv)
// a[3*ia1+idir1][3*ia2+idir2] =
// ( f[1][ia1][idir1][ia2][idir2] - f[0][ia1][idir1][ia2][idir2] ) / ( 2 h )
Context ctxt;
const int n = 3*nat;
DoubleMatrix a(ctxt,n,n);
valarray<double> a(n*n);
for ( int ia1 = 0; ia1 < nat; ia1++ )
{
......@@ -112,19 +112,25 @@ int main(int argc, char **argv)
double aij = ( f[1][ia1][idir1][ia2][idir2] -
f[0][ia1][idir1][ia2][idir2] ) / ( 2.0 * h );
aij = aij / sqrt(mass[ia1]*mass[ia2]);
a[j*a.n()+i] = aij;
a[j*n+i] = aij;
}
}
}
}
// cout << a;
valarray<double> w(a.m());
DoubleMatrix asave(a);
a.syev('u',w);
valarray<double> w(n);
valarray<double> asave(a);
char jobz = 'n';
char uplo = 'u';
int lwrk = 3*n;
valarray<double> wrk(lwrk);
int info;
dsyev_(&jobz,&uplo,&n,&a[0],&n,&w[0],&wrk[0],&lwrk,&info);
assert(info==0);
cout << " frequencies:" << endl;
for ( int i = 0; i < a.n(); i++ )
for ( int i = 0; i < n; i++ )
{
if ( w[i] < 0.0 )
cout << setw(8) << (int) (Ha2cm1 * sqrt(-w[i])) << " I" << endl;
......@@ -133,9 +139,11 @@ int main(int argc, char **argv)
}
a = asave;
a.syev('l',w);
dsyev_(&jobz,&uplo,&n,&a[0],&n,&w[0],&wrk[0],&lwrk,&info);
assert(info==0);
cout << " frequencies:" << endl;
for ( int i = 0; i < a.n(); i++ )
for ( int i = 0; i < n; i++ )
{
if ( w[i] < 0.0 )
cout << setw(8) << (int) (Ha2cm1 * sqrt(-w[i])) << " I" << endl;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment