Commit e894cae8 authored by eric pellegrini's avatar eric pellegrini
Browse files

Replaced Lapack dgesdd explicit call for SVD decomposition by numpy.linalg.svd

parent 736c2038
......@@ -11,18 +11,6 @@ import numpy
from MMTK import ParticleProperties, Visualization
dgesdd = None
try:
from numpy.linalg.lapack_lite import dgesdd
except ImportError:
pass
if dgesdd is None:
try:
# PyLAPACK
from lapack_dge import dgesdd
except ImportError: pass
#
# Class for a single mode
#
......@@ -229,64 +217,21 @@ class NormalModes(object):
ntotal = nexcluded + nmodes
natoms = len(basis[0])
sv = numpy.zeros((min(ntotal, 3*natoms),), numpy.float)
work_size = numpy.zeros((1,), numpy.float)
if nexcluded > 0:
self.basis = numpy.zeros((ntotal, 3*natoms), numpy.float)
for i in range(nexcluded):
self.basis[i] = numpy.ravel(excluded[i].array*self.weights)
min_n_m = min(3*natoms, nexcluded)
vt = numpy.zeros((nexcluded, min_n_m), numpy.float)
iwork = numpy.zeros((8*min_n_m,), int_type)
if 3*natoms >= nexcluded:
result = dgesdd('O', 3*natoms, nexcluded, self.basis,
3*natoms, sv, self.basis, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('O', 3*natoms, nexcluded, self.basis,
3*natoms, sv, self.basis, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
else:
u = numpy.zeros((min_n_m, 3*natoms), numpy.float)
result = dgesdd('S', 3*natoms, nexcluded, self.basis,
3*natoms, sv, u, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('S', 3*natoms, nexcluded, self.basis,
3*natoms, sv, u, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
self.basis[:min_n_m] = u
if result['info'] != 0:
raise ValueError('Lapack SVD: ' + `result['info']`)
_,sv,self.basis = numpy.linalg(self.basis,False)
svmax = numpy.maximum.reduce(sv)
nexcluded = numpy.add.reduce(numpy.greater(sv, 1.e-10*svmax))
ntotal = nexcluded + nmodes
for i in range(nmodes):
self.basis[i+nexcluded] = numpy.ravel(basis[i].array*self.weights)
min_n_m = min(3*natoms, ntotal)
vt = numpy.zeros((ntotal, min_n_m), numpy.float)
if 3*natoms >= ntotal:
result = dgesdd('O', 3*natoms, ntotal, self.basis, 3*natoms,
sv, self.basis, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
if int(work_size[0]) > work.shape[0]:
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('O', 3*natoms, ntotal, self.basis, 3*natoms,
sv, self.basis, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
else:
u = numpy.zeros((min_n_m, 3*natoms), numpy.float)
result = dgesdd('S', 3*natoms, ntotal, self.basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
if int(work_size[0]) > work.shape[0]:
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('S', 3*natoms, ntotal, self.basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
self.basis[:min_n_m] = u
if result['info'] != 0:
raise ValueError('Lapack SVD: ' + `result['info']`)
_,sv,self.basis = numpy.linalg(self.basis,False)
svmax = numpy.maximum.reduce(sv)
ntotal = numpy.add.reduce(numpy.greater(sv, 1.e-10*svmax))
nmodes = ntotal - nexcluded
......@@ -298,29 +243,9 @@ class NormalModes(object):
self.basis = numpy.array(map(lambda v: v.array, basis))
numpy.multiply(self.basis, self.weights, self.basis)
self.basis.shape = (nmodes, 3*natoms)
min_n_m = min(3*natoms, nmodes)
vt = numpy.zeros((nmodes, min_n_m), numpy.float)
iwork = numpy.zeros((8*min_n_m,), int_type)
if 3*natoms >= nmodes:
result = dgesdd('O', 3*natoms, nmodes, self.basis, 3*natoms,
sv, self.basis, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('O', 3*natoms, nmodes, self.basis, 3*natoms,
sv, self.basis, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
else:
u = numpy.zeros((min_n_m, 3*natoms), numpy.float)
result = dgesdd('S', 3*natoms, nmodes, self.basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work_size, -1, iwork, 0)
work = numpy.zeros((int(work_size[0]),), numpy.float)
result = dgesdd('S', 3*natoms, nmodes, self.basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
self.basis[:min_n_m] = u
if result['info'] != 0:
raise ValueError('Lapack SVD: ' + `result['info']`)
_,sv,self.basis = numpy.linalg(self.basis,False)
svmax = numpy.maximum.reduce(sv)
nmodes = numpy.add.reduce(numpy.greater(sv, 1.e-10*svmax))
ntotal = nmodes
......
......@@ -23,44 +23,6 @@ from MMTK import Utility, ParticleProperties
from Scientific.Geometry import Vector, ex, ey, ez
import copy
#
# Import LAPACK routines
#
dgesdd = None
try:
from numpy.linalg.lapack_lite import dgesdd, LapackError
except ImportError:
pass
if dgesdd is None:
try:
# PyLAPACK
from lapack_dge import dgesdd
except ImportError: pass
if dgesdd:
n = 1
array = numpy.zeros((n, n), numpy.float)
sv = numpy.zeros((n,), numpy.float)
u = numpy.zeros((n, n), numpy.float)
vt = numpy.zeros((n, n), numpy.float)
work = numpy.zeros((1,), numpy.float)
int_types = [numpy.int, numpy.int8, numpy.int16, numpy.int32]
try:
int_types.append(numpy.int64)
except AttributeError:
pass
for int_type in int_types:
iwork = numpy.zeros((1,), int_type)
try:
dgesdd('S', n, n, array, n, sv, u, n, vt, n, work, -1, iwork, 0)
break
except LapackError:
pass
del n, array, sv, u, vt, work, iwork, int_types
del array_package
#
# A set of particle vectors that define a subspace
#
......@@ -128,36 +90,15 @@ class Subspace(object):
nvectors = shape[0]
natoms = shape[1]
basis.shape = (nvectors, 3*natoms)
sv = numpy.zeros((min(nvectors, 3*natoms),), numpy.float)
min_n_m = min(3*natoms, nvectors)
vt = numpy.zeros((nvectors, min_n_m), numpy.float)
work = numpy.zeros((1,), numpy.float)
iwork = numpy.zeros((8*min_n_m,), int_type)
if 3*natoms >= nvectors:
result = dgesdd('O', 3*natoms, nvectors, basis, 3*natoms,
sv, basis, 3*natoms, vt, min_n_m,
work, -1, iwork, 0)
work = numpy.zeros((int(work[0]),), numpy.float)
result = dgesdd('O', 3*natoms, nvectors, basis, 3*natoms,
sv, basis, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
u = basis
else:
u = numpy.zeros((min_n_m, 3*natoms), numpy.float)
result = dgesdd('S', 3*natoms, nvectors, basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work, -1, iwork, 0)
work = numpy.zeros((int(work[0]),), numpy.float)
result = dgesdd('S', 3*natoms, nvectors, basis, 3*natoms,
sv, u, 3*natoms, vt, min_n_m,
work, work.shape[0], iwork, 0)
if result['info'] != 0:
raise ValueError('Lapack SVD: ' + `result['info']`)
_,sv,u = numpy.linalg.svd(basis,True)
svmax = numpy.maximum.reduce(sv)
nvectors = numpy.add.reduce(numpy.greater(sv, 1.e-10*svmax))
u = u[:nvectors]
u.shape = (nvectors, natoms, 3)
self._basis = ParticleVectorSet(self.universe, u)
return self._basis
def projectionOf(self, vector):
......
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