Write eigenvectors of collision matrices at temperatures.

    Depending on the choice of the solver, eigenvectors are sotred in
    either column-wise or row-wise.
170448
def write_unitary_matrix_to_hdf5(temperature,
                                 mesh,
                                 unitary_matrix=None,
                                 sigma=None,
                                 sigma_cutoff=None,
                                 solver=None,
                                 filename=None,
                                 verbose=False):
    

    suffix = _get_filename_suffix(mesh,
                                  sigma=sigma,
                                  sigma_cutoff=sigma_cutoff,
                                  filename=filename)
    hdf5_filename = "unitary" + suffix + ".hdf5"
    with h5py.File(hdf5_filename, 'w') as w:
        w.create_dataset('temperature', data=temperature)
        if unitary_matrix is not None:
            w.create_dataset('unitary_matrix', data=unitary_matrix)
        if solver is not None:
            w.create_dataset('solver', data=solver)

        if verbose:
            if len(temperature) > 1:
                text = "Unitary matrices "
            else:
                text = "Unitary matrix "
            if sigma is not None:
                text += "at sigma %s " % _del_zeros(sigma)
                if sigma_cutoff is not None:
                    text += "(%4.2f SD) " % sigma_cutoff
            if len(temperature) > 1:
                text += "were written into "
            else:
                text += "was written into "
            if sigma is not None:
                text += "\n"
            text += "\"%s\"." % hdf5_filename
            print(text)
324750
def eigendecompose(tensor, normalise=False):
    '''
    
    '''
    if normalise:
        tensor, tensor_norm = normalise_tensor(tensor)
    else:
        tensor_norm = 1.

    eigvals, eigvects = np.linalg.eigh(tensor, UPLO='U')

    isrt = np.argsort(eigvals)
    eigenvalues = eigvals[isrt] * tensor_norm
    eigenvectors = eigvects[:, isrt]
    return eigenvalues, eigenvectors
754252
def eigen(matrix):
    

    (vals, vecs) = np.linalg.eigh(matrix)
    ind = vals.argsort()[::-1]
    vals = vals[ind]
    vecs = vecs[:, ind]
    vals_ = vals.copy()
    vals_[vals_ < 0] = 0.
    cum_var_exp = np.cumsum(vals_ / vals_.sum())
    return Decomp(matrix.copy(), vals, vecs, cum_var_exp)
105926
def coarsen_matrix(Z, xlevel=0, ylevel=0, method='average'):
    

    # coarsen x
    if not ylevel:
        Z_coarsened = Z
    else:
        temp = []
        for z in Z: temp.append(coarsen_array(z, ylevel, method))
        Z_coarsened = _n.array(temp)

    # coarsen y
    if xlevel:
        Z_coarsened = Z_coarsened.transpose()
        temp = []
        for z in Z_coarsened: temp.append(coarsen_array(z, xlevel, method))
        Z_coarsened = _n.array(temp).transpose()

    return Z_coarsened



    # first coarsen the columns (if necessary)
    if ylevel:
        Z_ycoarsened = []
        for c in Z: Z_ycoarsened.append(coarsen_array(c, ylevel, method))
        Z_ycoarsened = _n.array(Z_ycoarsened)

    # now coarsen the rows
    if xlevel: return coarsen_array(Z_ycoarsened, xlevel, method)
    else:      return _n.array(Z_ycoarsened)
687389
private SparseDoubleVector getColumnVector(int column, boolean shouldLock) {
        int r = rows.get();
        if (shouldLock)
            lockColumn(column, r);
        // Ensure that the column data is up to date 
        while (lastVectorCacheUpdate.get() != modifications.get()) 
            updateVectorCache();
        int[] rowArr = colToRowsCache[column];
        SparseDoubleVector colVec = new SparseHashDoubleVector(r);
        for (int row : rowArr)
            colVec.set(row, matrixEntries.get(new Entry(row, column)));
        if (shouldLock)
            unlockColumn(column, r);
        return colVec;
    }
318007
private void setW( DMatrixRMaj W, DMatrixRMaj H , Vector3D_F64 a , Vector3D_F64 b ) {
		GeometryMath_F64.mult(H,b, tempV);
		setColumn(W,1, tempV);

		GeometryMath_F64.mult(H,a, tempV);
		setColumn(W,0, tempV);

		GeometryMath_F64.crossMatrix(tempV,Hv2);
		CommonOps_DDRM.mult(Hv2,H,tempM);
		GeometryMath_F64.mult(tempM,b, tempV);
		setColumn(W,2, tempV);
	}
45017
def _set_mode_kappa_Chaput(self, i_sigma, i_temp, weights):
        

        X = self._get_X(i_temp, weights, self._gv).ravel()
        num_ir_grid_points = len(self._ir_grid_points)
        num_band = self._primitive.get_number_of_atoms() * 3
        size = num_ir_grid_points * num_band * 3
        v = self._collision_matrix[i_sigma, i_temp].reshape(size, size)
        solver = _select_solver(self._pinv_solver)
        if solver in [1, 2, 4, 5]:
            v = v.T
        e = self._get_eigvals_pinv(i_sigma, i_temp)
        t = self._temperatures[i_temp]

        omega_inv = np.empty(v.shape, dtype='double', order='C')
        np.dot(v, (e * v).T, out=omega_inv)
        Y = np.dot(omega_inv, X)
        self._set_f_vectors(Y, num_ir_grid_points, weights)
        elems = ((0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1))
        for i, vxf in enumerate(elems):
            mat = self._get_I(vxf[0], vxf[1], num_ir_grid_points * num_band)
            self._mode_kappa[i_sigma, i_temp, :, :, i] = 0
            if mat is not None:
                np.dot(mat, omega_inv, out=mat)
                # vals = (X ** 2 * np.diag(mat)).reshape(-1, 3).sum(axis=1)
                # vals = vals.reshape(num_ir_grid_points, num_band)
                # self._mode_kappa[i_sigma, i_temp, :, :, i] = vals
                w = diagonalize_collision_matrix(mat,
                                                 pinv_solver=self._pinv_solver,
                                                 log_level=self._log_level)
                if solver in [1, 2, 4, 5]:
                    mat = mat.T
                spectra = np.dot(mat.T, X) ** 2 * w
                for s, eigvec in zip(spectra, mat.T):
                    vals = s * (eigvec ** 2).reshape(-1, 3).sum(axis=1)
                    vals = vals.reshape(num_ir_grid_points, num_band)
                    self._mode_kappa[i_sigma, i_temp, :, :, i] += vals

        factor = self._conversion_factor * Kb * t ** 2
        self._mode_kappa[i_sigma, i_temp] *= factor
581972
public void solve(DenseMatrix B) {
        int m = lu.nrows();
        int n = lu.ncols();
        int nrhs = B.ncols();

        if (B.nrows() != m)
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", lu.nrows(), lu.ncols(), B.nrows(), B.ncols()));

        if (isSingular()) {
            throw new RuntimeException("Matrix is singular.");
        }

        DenseMatrix X = Matrix.zeros(B.nrows(), B.ncols());

        // Copy right hand side with pivoting
        for (int j = 0; j < nrhs; j++) {
            for (int i = 0; i < m; i++) {
                X.set(i, j, B.get(piv[i], j));
            }
        }

        // Solve L*Y = B(piv,:)
        for (int k = 0; k < n; k++) {
            for (int i = k + 1; i < n; i++) {
                for (int j = 0; j < nrhs; j++) {
                    X.sub(i, j, X.get(k, j) * lu.get(i, k));
                }
            }
        }

        // Solve U*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                X.div(k, j, lu.get(k, k));
            }

            for (int i = 0; i < k; i++) {
                for (int j = 0; j < nrhs; j++) {
                    X.sub(i, j, X.get(k, j) * lu.get(i, k));
                }
            }
        }

        // Copy the result back to B.
        for (int j = 0; j < nrhs; j++) {
            for (int i = 0; i < m; i++) {
                B.set(i, j, X.get(i, j));
            }
        }
    }
1111590
private void gather(double[] z, SparseVector v, double taui, int d) {
        // Number of entries in the lower and upper part of the original matrix
        int nl = 0, nu = 0;
        for (VectorEntry e : v) {
            if (e.index() < d)
                nl++;
            else if (e.index() > d)
                nu++;
        }
        v.zero();

        // Entries in the L part of the vector
        lower.clear();
        for (int i = 0; i < d; ++i)
            if (Math.abs(z[i]) > taui)
                lower.add(new IntDoubleEntry(i, z[i]));

        // Entries in the U part of the vector
        upper.clear();
        for (int i = d + 1; i < z.length; ++i)
            if (Math.abs(z[i]) > taui)
                upper.add(new IntDoubleEntry(i, z[i]));

        // Sort in descending order
        Collections.sort(lower);
        Collections.sort(upper);

        // Always keep the diagonal
        v.set(d, z[d]);

        // Keep at most nl+p lower entries
        for (int i = 0; i < Math.min(nl + p, lower.size()); ++i) {
            IntDoubleEntry e = lower.get(i);
            v.set(e.index, e.value);
        }

        // Keep at most nu+p upper entries
        for (int i = 0; i < Math.min(nu + p, upper.size()); ++i) {
            IntDoubleEntry e = upper.get(i);
            v.set(e.index, e.value);
        }
    }
916924
def build_diagonals(self):
    

    ##########################################################
    # INCORPORATE BOUNDARY CONDITIONS INTO COEFFICIENT ARRAY #
    ##########################################################

    # Roll to keep the proper coefficients at the proper places in the
    # arrays: Python will naturally just do vertical shifts instead of 
    # diagonal shifts, so this takes into account the horizontal compoent 
    # to ensure that boundary values are at the right place.
    self.l2 = np.roll(self.l2, -2)
    self.l1 = np.roll(self.l1, -1)
    self.r1 = np.roll(self.r1, 1)
    self.r2 = np.roll(self.r2, 2)

    # Then assemble these rows: this is where the periodic boundary condition 
    # can matter.
    if self.coeff_matrix is not None:
      pass
    elif self.BC_E == 'Periodic' and self.BC_W == 'Periodic':
      # In this case, the boundary-condition-related stacking has already 
      # happened inside b.c.-handling function. This is because periodic
      # boundary conditions require extra diagonals to exist on the edges of 
      # the solution array
      pass
    else:
      self.diags = np.vstack((self.l2,self.l1,self.c0,self.r1,self.r2))
      self.offsets = np.array([-2,-1,0,1,2])

    # Everybody now (including periodic b.c. cases)
    self.coeff_matrix = spdiags(self.diags, self.offsets, self.nx, self.nx, format='csr')
49650
def _write_items(self, row_lookup, col_lookup, item):
        
        self.qc.write_items(row_lookup, col_lookup, item)
590725
public Matrix solve (Matrix B) {
      if (B.getRowDimension() != m) {
         throw new IllegalArgumentException("Matrix row dimensions must agree.");
      }
      if (!this.isNonsingular()) {
         throw new RuntimeException("Matrix is singular.");
      }

      // Copy right hand side with pivoting
      int nx = B.getColumnDimension();
      Matrix Xmat = B.getMatrix(piv,0,nx-1);
      double[][] X = Xmat.getArray();

      // Solve L*Y = B(piv,:)
      for (int k = 0; k < n; k++) {
         for (int i = k+1; i < n; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      // Solve U*X = Y;
      for (int k = n-1; k >= 0; k--) {
         for (int j = 0; j < nx; j++) {
            X[k][j] /= LU[k][k];
         }
         for (int i = 0; i < k; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      return Xmat;
   }
209090
public static void swapRows(Matrix matrix, long row1, long row2) {
		double temp = 0;
		long cols = matrix.getColumnCount();
		for (long col = 0; col < cols; col++) {
			temp = matrix.getAsDouble(row1, col);
			matrix.setAsDouble(matrix.getAsDouble(row2, col), row1, col);
			matrix.setAsDouble(temp, row2, col);
		}
	}
927880
public void solve(DenseMatrix B) {
        if (B.nrows() != L.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", L.nrows(), L.ncols(), B.nrows(), B.ncols()));
        }

        int n = B.nrows();
        int nrhs = B.ncols();

        // Solve L*Y = B;
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = 0; i < k; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(k, i));
                }
                B.div(k, j, L.get(k, k));
            }
        }

        // Solve L'*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = k + 1; i < n; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(i, k));
                }
                B.div(k, j, L.get(k, k));
            }
        }
    }
347553
function makeMatrix(columns, rows) {
  let matrix = [];
  for (let i = 0; i < columns; i++) {
    let column = [];
    for (let j = 0; j < rows; j++) {
      column.push(0);
    }
    matrix.push(column);
  }
  return matrix;
}
447524
public static ComplexFloatMatrix[] eigenvectors(FloatMatrix A) {
        A.assertSquare();
        // setting up result arrays
        FloatMatrix WR = new FloatMatrix(A.rows);
        FloatMatrix WI = WR.dup();
        FloatMatrix VR = new FloatMatrix(A.rows, A.rows);

        SimpleBlas.geev('N', 'V', A.dup(), WR, WI, dummyFloat, VR);

        // transferring the result
        ComplexFloatMatrix E = new ComplexFloatMatrix(WR, WI);
        ComplexFloatMatrix V = new ComplexFloatMatrix(A.rows, A.rows);
        //System.err.printf("VR = %s\n", VR.toString());
        for (int i = 0; i < A.rows; i++) {
            if (E.get(i).isReal()) {
                V.putColumn(i, new ComplexFloatMatrix(VR.getColumn(i)));
            } else {
                ComplexFloatMatrix v = new ComplexFloatMatrix(VR.getColumn(i), VR.getColumn(i + 1));
                V.putColumn(i, v);
                V.putColumn(i + 1, v.conji());
                i += 1;
            }
        }
        return new ComplexFloatMatrix[]{V, ComplexFloatMatrix.diag(E)};
    }
817778
public static void swapRows(double[][] matrix, int row1, int row2) {
		double[] temp = matrix[row1];
		matrix[row1] = matrix[row2];
		matrix[row2] = temp;
	}