From c688e23a9aeb85f072796281117ad6228dff7a48 Mon Sep 17 00:00:00 2001 From: Francois Gygi Date: Fri, 2 May 2014 22:00:34 +0000 Subject: [PATCH] use LAPACK function dsyev git-svn-id: http://qboxcode.org/svn/qb/trunk@1458 cba15fb0-1239-40c8-b417-11db7ca47a34 --- util/dynmat/D3vector.h | 168 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ util/dynmat/Makefile | 36 +++++++----------------------------- util/dynmat/dynmat.C | 30 +++++++++++++++++++----------- 3 files changed, 26 insertions(+), 208 deletions(-) delete mode 100644 util/dynmat/D3vector.h diff --git a/util/dynmat/D3vector.h b/util/dynmat/D3vector.h deleted file mode 100644 index 16cdc36..0000000 --- a/util/dynmat/D3vector.h +++ /dev/null @@ -1,168 +0,0 @@ -//////////////////////////////////////////////////////////////////////////////// -// -// 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 . -// -//////////////////////////////////////////////////////////////////////////////// -// -// 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 -#include -#include - -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 diff --git a/util/dynmat/Makefile b/util/dynmat/Makefile index 72bfe0e..45b0838 100644 --- a/util/dynmat/Makefile +++ b/util/dynmat/Makefile @@ -1,6 +1,6 @@ #------------------------------------------------------------------------------- # -# 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 . # #------------------------------------------------------------------------------- -# $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 diff --git a/util/dynmat/dynmat.C b/util/dynmat/dynmat.C index becdbed..f184c91 100644 --- a/util/dynmat/dynmat.C +++ b/util/dynmat/dynmat.C @@ -40,8 +40,9 @@ #include 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 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 w(a.m()); - DoubleMatrix asave(a); - a.syev('u',w); + valarray w(n); + valarray asave(a); + char jobz = 'n'; + char uplo = 'u'; + int lwrk = 3*n; + valarray 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; -- libgit2 0.26.0