cavis/libnd4j/include/helpers/cpu/biDiagonalUp.cpp

180 lines
5.6 KiB
C++
Raw Normal View History

2019-06-06 14:21:15 +02:00
/*******************************************************************************
* Copyright (c) 2015-2018 Skymind, Inc.
*
* This program and the accompanying materials are made available under the
* terms of the Apache License, Version 2.0 which is available at
* https://www.apache.org/licenses/LICENSE-2.0.
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
* License for the specific language governing permissions and limitations
* under the License.
*
* SPDX-License-Identifier: Apache-2.0
******************************************************************************/
//
// Created by Yurii Shyrma on 18.12.2017
//
#include <helpers/householder.h>
#include <helpers/biDiagonalUp.h>
#include <array/NDArrayFactory.h>
2019-06-06 14:21:15 +02:00
namespace sd {
2019-06-06 14:21:15 +02:00
namespace ops {
namespace helpers {
//////////////////////////////////////////////////////////////////////////
BiDiagonalUp::BiDiagonalUp(const NDArray& matrix): _HHmatrix(sd::NDArrayFactory::create(matrix.ordering(), {matrix.sizeAt(0), matrix.sizeAt(1)}, matrix.dataType(), matrix.getContext())),
_HHbidiag(sd::NDArrayFactory::create(matrix.ordering(), {matrix.sizeAt(1), matrix.sizeAt(1)}, matrix.dataType(), matrix.getContext())) {
2019-06-06 14:21:15 +02:00
// input validation
if(matrix.rankOf() != 2 || matrix.isScalar())
throw std::runtime_error("ops::helpers::biDiagonalizeUp constructor: input array must be 2D matrix !");
_HHmatrix.assign(&matrix);
_HHbidiag.assign(0.);
evalData();
}
template <typename T>
void BiDiagonalUp::_evalData() {
const auto rows = _HHmatrix.sizeAt(0);
const auto cols = _HHmatrix.sizeAt(1);
if(rows < cols)
throw std::runtime_error("ops::helpers::BiDiagonalizeUp::evalData method: this procedure is applicable only for input matrix with rows >= cols !");
NDArray* bottomRightCorner(nullptr), *column(nullptr), *row(nullptr);
T coeff, normX;
T _x, _y;
for(Nd4jLong i = 0; i < cols-1; ++i ) {
2019-06-06 14:21:15 +02:00
// evaluate Householder matrix nullifying columns
column = new NDArray(_HHmatrix({i,rows, i,i+1}, true));
_x = _HHmatrix.e<T>(i,i);
_y = _HHbidiag.e<T>(i,i);
Householder<T>::evalHHmatrixDataI(*column, _x, _y);
_HHmatrix.p<T>(i, i, _x);
_HHbidiag.p<T>(i, i, _y);
// multiply corresponding matrix block on householder matrix from the left: P * bottomRightCorner
bottomRightCorner = new NDArray(_HHmatrix({i,rows, i+1,cols}, true)); // {i, cols}
Householder<T>::mulLeft(*bottomRightCorner, _HHmatrix({i+1,rows, i,i+1}, true), _HHmatrix.e<T>(i,i));
delete bottomRightCorner;
delete column;
if(i == cols-2)
continue; // do not apply right multiplying at last iteration
// evaluate Householder matrix nullifying rows
row = new NDArray(_HHmatrix({i,i+1, i+1,cols}, true));
_x = _HHmatrix.e<T>(i,i+1);
_y = _HHbidiag.e<T>(i,i+1);
Householder<T>::evalHHmatrixDataI(*row, _x, _y);
_HHmatrix.p<T>(i, i+1, _x);
_HHbidiag.p<T>(i, i+1, _y);
// multiply corresponding matrix block on householder matrix from the right: bottomRightCorner * P
bottomRightCorner = new NDArray(_HHmatrix({i+1,rows, i+1,cols}, true)); // {i, rows}
Householder<T>::mulRight(*bottomRightCorner, _HHmatrix({i,i+1, i+2,cols}, true), _HHmatrix.e<T>(i,i+1));
delete bottomRightCorner;
delete row;
}
row = new NDArray(_HHmatrix({cols-2,cols-1, cols-1,cols}, true));
_x = _HHmatrix.e<T>(cols-2,cols-1);
_y = _HHbidiag.e<T>(cols-2,cols-1);
Householder<T>::evalHHmatrixDataI(*row, _x, _y);
_HHmatrix.p<T>(cols-2,cols-1, _x);
_HHbidiag.p<T>(cols-2,cols-1, _y);
delete row;
column = new NDArray(_HHmatrix({cols-1,rows, cols-1,cols}, true));
_x = _HHmatrix.e<T>(cols-1,cols-1);
_y = _HHbidiag.e<T>(cols-1,cols-1);
Householder<T>::evalHHmatrixDataI(*column, _x, _y);
_HHmatrix.p<T>(cols-1, cols-1, _x);
_HHbidiag.p<T>(cols-1, cols-1, _y);
delete column;
}
//////////////////////////////////////////////////////////////////////////
void BiDiagonalUp::evalData() {
auto xType = _HHmatrix.dataType();
BUILD_SINGLE_SELECTOR(xType, _evalData, ();, FLOAT_TYPES);
}
//////////////////////////////////////////////////////////////////////////
template <typename T>
HHsequence BiDiagonalUp::makeHHsequence_(const char type) const {
if(type == 'u') {
const int diagSize = _HHbidiag.sizeAt(0);
auto colOfCoeffs = NDArrayFactory::create(_HHmatrix.ordering(), {diagSize, 1}, _HHmatrix.dataType(), _HHmatrix.getContext());
for(int i = 0; i < diagSize; ++i)
colOfCoeffs.p(i, _HHmatrix.e<T>(i,i));
return HHsequence(_HHmatrix, colOfCoeffs, type);
}
else {
const int diagUpSize = _HHbidiag.sizeAt(0) - 1;
NDArray colOfCoeffs = NDArrayFactory::create(_HHmatrix.ordering(), {diagUpSize, 1}, _HHmatrix.dataType(), _HHmatrix.getContext());
for(int i = 0; i < diagUpSize; ++i)
colOfCoeffs.p(i, _HHmatrix.e<T>(i,i+1));
HHsequence result(_HHmatrix, colOfCoeffs, type);
result._diagSize = diagUpSize;
result._shift = 1;
return result;
}
}
HHsequence BiDiagonalUp::makeHHsequence(const char type) const {
auto xType = _HHmatrix.dataType();
BUILD_SINGLE_SELECTOR(xType, return makeHHsequence_, (type);, FLOAT_TYPES);
}
BUILD_SINGLE_TEMPLATE(template void BiDiagonalUp::_evalData, (), FLOAT_TYPES);
BUILD_SINGLE_TEMPLATE(template HHsequence BiDiagonalUp::makeHHsequence_, (const char type) const, FLOAT_TYPES);
}
}
}