GPAW’s Matrix object

A simple example that we can run with MPI on 4 cores:

from gpaw.matrix import Matrix
from gpaw.mpi import world
a = Matrix(5, 5, dist=(world, 2, 2, 2))
a.array[:] = world.rank
print(world.rank, a.array.shape)

Here, we have created a 5x5 Matrix of floats distributed on a 2x2 BLACS gris with a blocksize 2 and we then print the shapes of the ndarrays, which looks like this (in random order):

1 (2, 3)
2 (3, 2)
3 (2, 2)
0 (3, 3)

Let’s create a new matrix b and redistribute from a to b:

b = a.new(dist=(None, 1, 1, None))
a.redist(b)
if world.rank == 0:
    print(b.array)

This will output:

[[ 0.  0.  2.  2.  0.]
 [ 0.  0.  2.  2.  0.]
 [ 1.  1.  3.  3.  1.]
 [ 1.  1.  3.  3.  1.]
 [ 0.  0.  2.  2.  0.]]

Matrix-matrix multiplication works like this:

from gpaw.matrix import matrix_matrix_multiply as mmm
c = mmm(1.0, a, 'N', a, 'T')
mmm(1.0, a, 'N', a, 'T', 1.0, c, symmetric=True)

gpaw.matrix.matrix_matrix_multiply(alpha, a, opa, b, opb, beta=0.0, c=None, symmetric=False)

BLAS-style matrix-matrix multiplication.

Will use dgemm/zgemm/dsyrk/zherk/dsyr2k/zher2k as apropriate or the equivalent PBLAS functions for distributed matrices.

The coefficients alpha and beta are of type float. Matrices a, b and c must have same type (float or complex). The strings apa and opb must be ‘N’, ‘T’, or ‘C’ . For opa=’N’ and opb=’N’, the operation performed is equivalent to:

c.array[:] =  alpha * np.dot(a.array, b.array) + beta * c.array

Replace a.array with a.array.T or a.array.T.conj() for opa=’T’ and ‘C’ resprctively (similarly for opb).

Use symmetric=True if the result matrix is symmetric/hermetian (only lower half of c will be evaluated).

class gpaw.matrix.Matrix(M, N, dtype=None, data=None, dist=None)

Matrix object.

M: int

Rows.

N: int

Columns.

dtype: type

Data type (float or complex).

dist: tuple or None

BLACS distribution given as (communicator, rows, colums, blocksize) tuple. Default is None meaning no distribution.

data: ndarray or None.

Numpy ndarray to use for starage. By default, a new ndarray will be allocated.

complex_conjugate()

Inplace complex conjugation.

eigh(cc=False, scalapack=(None, 1, 1, None))

Calculate eigenvectors and eigenvalues.

Matrix must be symmetric/hermitian and stored in lower half.

cc: bool

Complex conjugate matrix before finding eigenvalues.

scalapack: tuple

BLACS distribution for ScaLapack to use. Default is to do serial diagonalization.

invcholesky()

Inverse of Cholesky decomposition.

Only the lower part is used.

multiply(alpha, opa, b, opb, beta=0.0, out=None, symmetric=False)

BLAS-style Matrix-matrix multiplication.

See matrix_matrix_multipliction() for details.

new(dist='inherit')

Create new matrix of same shape and dtype.

Default is to use same BLACS distribution. Use dist to use another distribution.

redist(other)

Redistribute to other BLACS layout.