Clearly this transformation, in general, mixes data between all of the nodes. Rather than storing an entire wavefunction we transform just the parts needed for one node at a time.
! construct the data for node `node' do node=1,bnd_nodes transform local data to give contribution to `node' reduce over bnd-group end do
This is simplest if we encapsulate all such operations into a single subroutine, overloaded for the various wavefunction types. In fact such a routine already exists, wave_rotate.
There are also occasions when we wish to transform a wavefunction by a matrix that is not square, so that the resulting number of bands is not the same as the initial number of bands. This occurs when, for example, we orthogonalise a slice to an existing wavefunction; we copy the wavefunction to a temporary slice, transform the data by multiplying by the overlap matrix, which changes the number of bands to the same as the slice, and then subtract it from slice. This is not a `rotation', and goes against the specification for wave_rotate. For this reason we coded a new, more general subroutine wave_transform which allows any linear transformation of the bands.
The orthogonalisation subroutines can now be coded compactly as:
! orthogonalise slice to all bands of wvfn call wave_dot_all(wvfn,slice,overlap) call wave_copy(wvfn,tmp_slice) call wave_transform(tmp_slice,overlap) ! Subtract overlapping components call wave_add(tmp_slice,slice,c1=-cmplx_1)