You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mxnet.apache.org by GitBox <gi...@apache.org> on 2019/12/28 03:13:07 UTC
[GitHub] [incubator-mxnet] sjtuWangDing opened a new pull request #17188:
[Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
sjtuWangDing opened a new pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188
## Description ##
Add linalg.eig/eigh/eigvals/eigvalsh op
## Checklist ##
### Essentials ###
Please feel free to remove inapplicable items for your PR.
- [ ] The PR title starts with [MXNET-$JIRA_ID], where $JIRA_ID refers to the relevant [JIRA issue](https://issues.apache.org/jira/projects/MXNET/issues) created (except PRs with tiny changes)
- [ ] Changes are complete (i.e. I finished coding on this PR)
- [ ] All changes have test coverage:
- Unit tests are added for small changes to verify correctness (e.g. adding a new operator)
- Nightly tests are added for complicated/long-running ones (e.g. changing distributed kvstore)
- Build tests will be added for build configuration changes (e.g. adding a new build option with NCCL)
- [ ] Code is well-documented:
- For user-facing API changes, API doc string has been updated.
- For new C++ functions in header files, their functionalities and arguments are documented.
- For new examples, README.md is added to explain the what the example does, the source of the dataset, expected performance on test set and reference to the original paper if applicable
- Check the API doc at https://mxnet-ci-doc.s3-accelerate.dualstack.amazonaws.com/PR-$PR_ID/$BUILD_ID/index.html
- [ ] To the best of my knowledge, examples are either not affected by this change, or have been fixed to be compatible with this change
### Changes ###
- [ ] Feature1, tests, (and when applicable, API doc)
- [ ] Feature2, tests, (and when applicable, API doc)
## Comments ##
- If this change is a backward incompatible change, why must this change be made.
- Interesting edge cases to note here
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 merged pull request #17188: [Numpy] Add
linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 merged pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364467821
##########
File path: src/operator/numpy/linalg/np_eigvals.cc
##########
@@ -0,0 +1,124 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals.cc
+ * \brief CPU implementation placeholder of Eigvals Operator
+ */
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig.
+bool EigvalsOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ CHECK_GE(eig_ndim, 1)
+ << "Outputs W must be at least one-dimensional";
+ std::vector<int> a_shape_vec(eig_ndim + 1);
+ for (int i = 0; i < eig_ndim; ++i) {
+ a_shape_vec[i] = eig_shape[i];
+ }
+ a_shape_vec[eig_ndim] = eig_shape[eig_ndim - 1];
+ mxnet::TShape a_shape(a_shape_vec.begin(), a_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, a_shape);
+ }
+ }
+ return shape_is_known(*in_attrs) && shape_is_known(*out_attrs);
+}
+
+inline bool EigvalsOpType(const nnvm::NodeAttrs& attrs,
+ std::vector<int>* in_attrs,
Review comment:
alignment
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364998737
##########
File path: python/mxnet/symbol/numpy/linalg.py
##########
@@ -496,3 +497,265 @@ def tensorsolve(a, b, axes=None):
True
"""
return _npi.tensorsolve(a, b, axes)
+
+
+def eigvals(a):
+ r"""
+ Compute the eigenvalues of a general matrix.
+
+ Main difference between `eigvals` and `eig`: the eigenvectors aren't
+ returned.
+
+ Parameters
+ ----------
+ a : (..., M, M) ndarray
+ A real-valued matrix whose eigenvalues will be computed.
+
+ Returns
+ -------
+ w : (..., M,) ndarray
+ The eigenvalues, each repeated according to its multiplicity.
+ They are not necessarily ordered.
+
+ Raises
+ ------
+ MXNetError
+ If the eigenvalue computation does not converge.
+
+ See Also
+ --------
+ eig : eigenvalues and right eigenvectors of general arrays
+ eigh : eigenvalues and eigenvectors of a real symmetric array.
+ eigvalsh : eigenvalues of a real symmetric.
+
+ Notes
+ -----
+ Broadcasting rules apply, see the `numpy.linalg` documentation for
+ details.
+
+ This is implemented using the ``_geev`` LAPACK routines which compute
+ the eigenvalues and eigenvectors of general square arrays.
+
+ This function differs from the original `numpy.linalg.eigvals
+ <https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigvals.html>`_ in
+ the following way(s):
+ - Does not support complex input and output.
+
+ Examples
+ --------
+ Illustration, using the fact that the eigenvalues of a diagonal matrix
+ are its diagonal elements, that multiplying a matrix on the left
+ by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
+ of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
+ if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
+ ``A``:
+ >>> from numpy import linalg as LA
+ >>> x = np.random.random()
+ >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
+ >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
+ (1.0, 1.0, 0.0)
+
+ Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
+ >>> D = np.diag((-1,1))
Review comment:
No need for examples for symbolic APIs.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] sjtuWangDing commented on a change in pull
request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
sjtuWangDing commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r366190317
##########
File path: src/operator/numpy/linalg/np_eigvals.cc
##########
@@ -0,0 +1,124 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals.cc
+ * \brief CPU implementation placeholder of Eigvals Operator
+ */
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig.
+bool EigvalsOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ CHECK_GE(eig_ndim, 1)
+ << "Outputs W must be at least one-dimensional";
+ std::vector<int> a_shape_vec(eig_ndim + 1);
+ for (int i = 0; i < eig_ndim; ++i) {
+ a_shape_vec[i] = eig_shape[i];
+ }
+ a_shape_vec[eig_ndim] = eig_shape[eig_ndim - 1];
+ mxnet::TShape a_shape(a_shape_vec.begin(), a_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, a_shape);
+ }
+ }
+ return shape_is_known(*in_attrs) && shape_is_known(*out_attrs);
+}
+
+inline bool EigvalsOpType(const nnvm::NodeAttrs& attrs,
+ std::vector<int>* in_attrs,
+ std::vector<int>* out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ int a_type = in_attrs->at(0);
+ // unsupport float16
Review comment:
These 2 ops support integer inputs in numpy.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r365001641
##########
File path: src/operator/numpy/linalg/np_eigvals.cc
##########
@@ -0,0 +1,124 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals.cc
+ * \brief CPU implementation placeholder of Eigvals Operator
+ */
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig.
+bool EigvalsOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ CHECK_GE(eig_ndim, 1)
+ << "Outputs W must be at least one-dimensional";
+ std::vector<int> a_shape_vec(eig_ndim + 1);
+ for (int i = 0; i < eig_ndim; ++i) {
+ a_shape_vec[i] = eig_shape[i];
+ }
+ a_shape_vec[eig_ndim] = eig_shape[eig_ndim - 1];
+ mxnet::TShape a_shape(a_shape_vec.begin(), a_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, a_shape);
+ }
+ }
+ return shape_is_known(*in_attrs) && shape_is_known(*out_attrs);
+}
+
+inline bool EigvalsOpType(const nnvm::NodeAttrs& attrs,
+ std::vector<int>* in_attrs,
+ std::vector<int>* out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ int a_type = in_attrs->at(0);
+ // unsupport float16
Review comment:
Does the 2 ops support integer inputs?
If not, then I would suggest the following code structure:
```c++
if (common::is_float(a_type)) {
// raise error for fp16
TYPE_ASSIGN_CHECK(*out_attrs, 0, in_attrs->at(0));
} else {
// raise error for integer
}
```
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] reminisce commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
reminisce commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r362083569
##########
File path: tests/python/unittest/test_numpy_op.py
##########
@@ -4291,6 +4291,377 @@ def newInvertibleMatrix_2D(shape, max_cond=4):
check_tensorsolve(mx_out, a.asnumpy(), b.asnumpy(), axes)
+@with_seed()
+@use_np
+def test_np_linalg_eigvals():
+ class TestEigvals(HybridBlock):
+ def __init__(self):
+ super(TestEigvals, self).__init__()
+
+ def hybrid_forward(self, F, a):
+ return F.np.linalg.eigvals(a)
+
+ def check_eigvals(x, a_np):
+ try:
+ x_expected = _np.linalg.eigvals(a_np)
+ except Exception as e:
+ print("a:", a_np)
+ print("a shape:", a_np.shape)
+ print(e)
+ else:
+ assert x.shape == x_expected.shape
+ if 0 not in x.shape:
+ n = int(_np.prod(x.shape[:-1])) if len(shape) > 1 else 1
+ x = x.reshape(n, -1)
+ x_expected = x_expected.reshape(n, -1)
+ for i in range(n):
+ x1 = _np.sort(x[i].asnumpy())
+ x2 = _np.sort(x_expected[i])
+ assert_almost_equal(x1, x2, rtol=rtol, atol=atol)
+
+ def newMatrixWithRealEigvals_2D(n):
+ shape = (n, n)
+ q = _np.ones(shape)
+ while 1:
+ D = _np.diag(_np.random.uniform(-1.0, 1.0, shape[-1]))
+ I = _np.eye(shape[-1]).reshape(shape)
+ v = _np.random.uniform(-1., 1., shape[-1]).reshape(shape[:-1] + (1,))
+ v = v / _np.linalg.norm(v, axis=-2, keepdims=True)
+ v_T = _np.swapaxes(v, -1, -2)
+ U = I - 2 * _np.matmul(v, v_T)
+ q = _np.matmul(U, D)
+ if (_np.linalg.cond(q, 2) < 3):
+ break
+ D = _np.diag(_np.random.uniform(-10.0, 10.0, n))
+ q_inv = _np.linalg.inv(q)
+ return _np.matmul(_np.matmul(q_inv, D), q)
+
+ def newMatrixWithRealEigvals_nD(shape):
+ n = int(_np.prod(shape[:-2])) if len(shape) > 2 else 1
+ return _np.array([newMatrixWithRealEigvals_2D(shape[-1]) for i in range(n)]).reshape(shape)
+
+ shapes = [
+ (0, 0),
+ (1, 1),
+ (3, 3),
+ (5, 5),
+ (1, 0, 0),
+ (0, 4, 4),
+ (1, 4, 4),
+ (2, 4, 4),
+ (5, 5, 5),
+ (1, 1, 4, 4),
+ (2, 3, 4, 4)
+ ]
+ dtypes = ['float32', 'float64']
+ for hybridize in [True, False]:
+ for shape, dtype in itertools.product(shapes, dtypes):
+ rtol = 1e-2 if dtype == 'float32' else 1e-3
+ atol = 1e-4 if dtype == 'float32' else 1e-5
+ test_eigvals = TestEigvals()
+ if hybridize:
+ test_eigvals.hybridize()
+ if 0 in shape:
+ a_np = _np.ones(shape)
+ else:
+ a_np = newMatrixWithRealEigvals_nD(shape)
+ a = np.array(a_np, dtype=dtype)
+ # check eig validity
+ w = test_eigvals(a)
+ check_eigvals(w, a.asnumpy())
+
+
+@with_seed()
+@use_np
+def test_np_linalg_eigvalsh():
+ class TestEigvalsh(HybridBlock):
+ def __init__(self, UPLO):
+ super(TestEigvalsh, self).__init__()
+ self._UPLO = UPLO
+
+ def hybrid_forward(self, F, a):
+ return F.np.linalg.eigvalsh(a, UPLO=self._UPLO)
+
+ def check_eigvalsh(w, a_np, UPLO):
+ try:
+ w_expected = _np.linalg.eigvalsh(a_np, UPLO)
+ except Exception as e:
+ print("a:", a_np)
+ print("a shape:", a_np.shape)
+ print(e)
+ else:
+ assert w.shape == w_expected.shape
+ assert_almost_equal(w, w_expected, rtol=rtol, atol=atol)
+
+ def newOrthonormalMatrix_2D(n):
Review comment:
Please add util functions to `mxnet.test_utils`. Also, please adopt python function's naming convention.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364468214
##########
File path: src/operator/numpy/linalg/np_eigvals-inl.h
##########
@@ -0,0 +1,432 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals-inl.h
+ * \brief Placeholder for eigvals
+ */
+#ifndef MXNET_OPERATOR_NUMPY_LINALG_NP_EIGVALS_INL_H_
+#define MXNET_OPERATOR_NUMPY_LINALG_NP_EIGVALS_INL_H_
+
+#include <mxnet/operator_util.h>
+#include <vector>
+#include "../../operator_common.h"
+#include "../../mshadow_op.h"
+#include "../../tensor/la_op.h"
+#include "../../tensor/la_op-inl.h"
+#include "./np_solve-inl.h"
+
+namespace mxnet {
+namespace op {
+
+using namespace mshadow;
+
+template<int req>
+struct eigvals_assign_helper {
+ template<typename DType>
+ MSHADOW_XINLINE static void Map(int i, const DType *in_data, DType *out_data) {
+ KERNEL_ASSIGN(out_data[i], req, in_data[i]);
+ }
+};
+
+struct eigh_eigvalsh_helper {
+ template<typename InDType, typename OutDType>
+ MSHADOW_XINLINE static void Map(int i, const InDType *in_data, OutDType *out_data,
+ const int nrow, const int ld, const int step, bool USE_UP) {
+ int idx = i / step, row = (i % step) / ld, col = (i % step) % ld;
+ if ((USE_UP && row > col) || (!USE_UP && row < col)) {
+ out_data[idx * step + col + row * ld] =
+ static_cast<OutDType>(in_data[idx * step + row + col * ld]);
+ } else {
+ out_data[idx * step + col + row * ld] =
+ static_cast<OutDType>(in_data[idx * step + col + row * ld]);
+ }
+ }
+};
+
+template<typename xpu, typename DType>
+void linalg_geev(char jobvl,
+ char jobvr,
+ const Tensor<xpu, 2, DType>& a,
+ const Tensor<xpu, 1, DType>& wr,
+ const Tensor<xpu, 1, DType>& wi,
+ const Tensor<xpu, 2, DType>& vl,
+ const Tensor<xpu, 2, DType>& vr,
+ const Tensor<xpu, 1, DType>& work_array);
+
+#define LINALG_CPU_EIG(fname, DType) \
+template<> inline \
+void linalg_geev<cpu, DType>(char jobvl, \
+ char jobvr, \
+ const Tensor<cpu, 2, DType>& a, \
+ const Tensor<cpu, 1, DType>& wr, \
+ const Tensor<cpu, 1, DType>& wi, \
+ const Tensor<cpu, 2, DType>& vl, \
+ const Tensor<cpu, 2, DType>& vr, \
+ const Tensor<cpu, 1, DType>& work_array) { \
+ const int n = a.size(1), lda = a.size(0); \
+ const int lwork = work_array.shape_.Size(); \
+ const int ldvl = vl.size(0), ldvr = vr.size(0); \
+ int res(MXNET_LAPACK_##fname(MXNET_LAPACK_COL_MAJOR, jobvl, jobvr, \
+ n, a.dptr_, lda, \
+ wr.dptr_, wi.dptr_, \
+ vl.dptr_, ldvl, \
+ vr.dptr_, ldvr, \
+ work_array.dptr_, lwork)); \
+ CHECK_LE(res, 0) << #fname << "the QR algorithm failed to compute all the" \
+ << "eigenvalues, and no eigenvectors have been computed; elements " \
+ << res + 1 << ":N" << " of WR and WI contain eigenvalues which have converged"; \
+ CHECK_GE(res, 0) << #fname << ": the " << -res \
+ << "-th argument had an illegal value"; \
+}
+LINALG_CPU_EIG(sgeev, float)
+LINALG_CPU_EIG(dgeev, double)
+
+#ifdef __CUDACC__
+
+#define LINALG_GPU_EIG(fname, DType) \
+template<> inline \
+void linalg_geev<gpu, DType>(char jobvl, \
+ char jobvr, \
+ const Tensor<gpu, 2, DType>& a, \
+ const Tensor<gpu, 1, DType>& wr, \
+ const Tensor<gpu, 1, DType>& wi, \
+ const Tensor<gpu, 2, DType>& vl, \
+ const Tensor<gpu, 2, DType>& vr, \
+ const Tensor<gpu, 1, DType>& work_array) { \
+ LOG(FATAL) << "Lapack _geev routines in gpu is unsupported"; \
+}
+LINALG_GPU_EIG(sgeev, float)
+LINALG_GPU_EIG(dgeev, double)
+
+#endif // __CUDACC__
+
+struct eig_eigvals {
+ template<typename xpu, typename DType>
+ static void op(char jobvl,
+ char jobvr,
+ const Tensor<xpu, 3, DType>& a,
+ const Tensor<xpu, 2, DType>& wr,
+ const Tensor<xpu, 2, DType>& wi,
+ const Tensor<xpu, 3, DType>& vl,
+ const Tensor<xpu, 3, DType>& vr,
+ const Tensor<xpu, 1, DType>& work_array) {
+ const mxnet::TShape& a_shape = a.shape_;
+ const int a_ndim = a_shape.ndim();
+ if (jobvl == 'N' && jobvr == 'N') {
+ CHECK_GE(work_array.shape_.Size(), 3 * a.shape_[a_ndim - 1])
+ << "The dimension of the array WORK in LAPACKE_#GEEV should >= max(1,3*N).";
+ } else {
+ CHECK_GE(work_array.shape_.Size(), 4 * a.shape_[a_ndim - 1])
+ << "If JOBVL = 'V' or JOBVR = 'V', "
+ << "the dimension of the array WORK in LAPACKE_#GEEV should >= 4*N.";
+ }
+ for (int i = 0; i < a_shape[0]; ++i) {
+ if (jobvl == 'N' && jobvr == 'N') {
+ linalg_geev(jobvl, jobvr, a[i], wr[i], wi[i], vl[0], vr[0], work_array);
+ } else if (jobvl == 'N' && jobvr == 'V') {
+ linalg_geev(jobvl, jobvr, a[i], wr[i], wi[i], vl[0], vr[i], work_array);
+ } else if (jobvl == 'V' && jobvr == 'N') {
+ linalg_geev(jobvl, jobvr, a[i], wr[i], wi[i], vl[i], vr[0], work_array);
+ } else {
+ linalg_geev(jobvl, jobvr, a[i], wr[i], wi[i], vl[i], vr[i], work_array);
+ }
+ }
+ }
+};
+
+// Calculates workspace size of eigvals forward op.
+// The dimension of the array WORK in LAPACKE_#GEEV should >= max(1,3*N), and
+// if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.
+// For good performance, LWORK must generally be larger.
+template<typename xpu>
+size_t EigvalsForwardWorkspaceSize(const TBlob& a,
+ const TBlob& w,
+ const std::vector<OpReqType>& req) {
+ if (kNullOp == req[0]) { return 0U; }
+
+ // Zero-size input, no need to launch kernel
+ if (0U == a.Size()) { return 0U; }
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, DType, {
+ size_t work_space_size = 0;
+ size_t n = a.size(a.ndim() - 1);
+ work_space_size += a.Size(); // For matrix.
+ work_space_size += 2 * w.Size(); // For eigenvalues' real and image component.
+ work_space_size += 2 * n * n; // For left and right eigenvectors temp memory
+ work_space_size += 3 * n; // For workspace size in LAPACKE_#GEEV.
+ work_space_size *= sizeof(DType);
+ return work_space_size;
+ });
+ LOG(FATAL) << "InternalError: cannot reach here";
+ return 0U;
+}
+
+template<typename xpu>
+void EigvalsOpForwardImpl(const TBlob& a,
+ const TBlob& w,
+ const std::vector<OpReqType>& req,
+ std::vector<char> *workspace,
+ mshadow::Stream<xpu> *s) {
+ if (kNullOp == req[0]) { return; }
+ const mxnet::TShape& a_shape = a.shape_;
+ const mxnet::TShape& w_shape = w.shape_;
+ const int a_ndim = a_shape.ndim();
+
+ // Zero-size output, no need to launch kernel
+ if (0U == a.Size()) { return; }
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, DType, {
+ const int N = a.size(a_ndim - 1);
+ DType *a_ptr =
+ reinterpret_cast<DType*>(workspace->data());
+ DType *wr_ptr =
+ reinterpret_cast<DType*>(workspace->data() + a.Size() * sizeof(DType));
+ DType *wi_ptr =
+ reinterpret_cast<DType*>(workspace->data() + (w.Size() + a.Size()) * sizeof(DType));
+ DType *vl_ptr =
+ reinterpret_cast<DType*>(workspace->data() + (2 * w.Size() + a.Size()) * sizeof(DType));
+ DType *vr_ptr =
+ reinterpret_cast<DType*>(
+ workspace->data() + (N * N + 2 * w.Size() + a.Size()) * sizeof(DType));
+ DType *work_ptr =
+ reinterpret_cast<DType*>(
+ workspace->data() + (2 * (N * N + w.Size()) + a.Size()) * sizeof(DType));
+ MSHADOW_TYPE_SWITCH(a.type_flag_, AType, {
+ // Cast type and transpose.
+ mxnet_op::Kernel<SolveTypeTransposeHelper, xpu>::Launch(
+ s, a.Size(), a.dptr<AType>(), a_ptr, N, N, N * N);
+ });
+ s->Wait();
+ char jobvl = 'N', jobvr = 'N';
+ mxnet::TBlob a_trans_data(a_ptr, a_shape, a.dev_mask(), a.dev_id());
+ mxnet::TBlob wr_data(wr_ptr, w_shape, w.dev_mask(), w.dev_id());
+ mxnet::TBlob wi_data(wi_ptr, w_shape, w.dev_mask(), w.dev_id());
+ mxnet::TBlob vl_data(vl_ptr, Shape3(1, N, N), w.dev_mask(), w.dev_id());
+ mxnet::TBlob vr_data(vr_ptr, Shape3(1, N, N), w.dev_mask(), w.dev_id());
+ mxnet::TBlob work_data(work_ptr, Shape1(3 * N), a.dev_mask(), a.dev_id());
+ eig_eigvals::op(jobvl, jobvr,
+ a_trans_data.FlatToKD<xpu, 3, DType>(s),
+ wr_data.FlatToKD<xpu, 2, DType>(s),
+ wi_data.FlatToKD<xpu, 2, DType>(s),
+ vl_data.get<xpu, 3, DType>(s),
+ vr_data.get<xpu, 3, DType>(s),
+ work_data.get<xpu, 1, DType>(s));
+ for (size_t i = 0; i < wi_data.Size(); ++i) {
+ CHECK_LE(fabs(wi_ptr[i]), 1e-15)
+ << "Complex eigvals is unsupported in linalg temporary.";
+ }
+ MXNET_ASSIGN_REQ_SWITCH(req[0], req_type, {
+ mxnet_op::Kernel<eigvals_assign_helper<req_type>, xpu>::Launch(
+ s, w.Size(), wr_ptr, w.dptr<DType>());
+ });
+ });
+}
+
+template<typename AType, typename WType>
+void GpuCallbackCpuImpl(const TBlob& a,
+ const TBlob& w,
+ AType* a_cp_ptr,
+ WType* w_cp_ptr,
+ std::vector<char>* workspace,
+ const OpContext& ctx,
+ const std::vector<OpReqType>& req) {
+#if MXNET_USE_CUDA
+ mshadow::Stream<gpu> *s = ctx.get_stream<gpu>();
+ cudaStream_t stream = Stream<gpu>::GetStream(s);
+ CUDA_CALL(cudaMemcpyAsync(a_cp_ptr, a.dptr<AType>(), sizeof(AType) * a.Size(),
+ cudaMemcpyDeviceToHost, stream));
+ CUDA_CALL(cudaMemcpyAsync(w_cp_ptr, w.dptr<WType>(), sizeof(WType) * w.Size(),
+ cudaMemcpyDeviceToHost, stream));
+ CUDA_CALL(cudaStreamSynchronize(stream));
+ mxnet::TBlob a_data(a_cp_ptr, a.shape_, cpu::kDevMask);
+ mxnet::TBlob w_data(w_cp_ptr, w.shape_, cpu::kDevMask);
+ // Op forward implement on cpu.
+ EigvalsOpForwardImpl<cpu>(a_data, w_data, req, workspace, ctx.get_stream<cpu>());
+ // Copy back to gpu.
+ CUDA_CALL(cudaMemcpyAsync(w.dptr<WType>(), w_cp_ptr, sizeof(WType) * w.Size(),
+ cudaMemcpyHostToDevice, stream));
+ CUDA_CALL(cudaStreamSynchronize(stream));
+#else
+ LOG(FATAL) << "Please build with USE_CUDA=1 to enable GPU";
+#endif // MXNET_USE_CUDA
+}
+
+template<typename xpu>
+void EigvalsOpForward(const nnvm::NodeAttrs& attrs,
+ const OpContext& ctx,
+ const std::vector<TBlob>& inputs,
+ const std::vector<OpReqType>& req,
+ const std::vector<TBlob>& outputs) {
+ CHECK_EQ(inputs.size(), 1U);
+ CHECK_EQ(outputs.size(), 1U);
+ CHECK_EQ(req.size(), 1U);
+ const TBlob& a = inputs[0];
+ const TBlob& w = outputs[0];
+
+ // Calculate workspace size.
+ size_t workspace_size = EigvalsForwardWorkspaceSize<cpu>(a, w, req);
+ std::vector<char> workspace(workspace_size, 0);
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, WType, {
+ MSHADOW_SGL_DBL_TYPE_SWITCH(a.type_flag_, AType, {
+ if (xpu::kDevCPU) {
+ // Op forward implement.
+ EigvalsOpForwardImpl<cpu>(a, w, req, &workspace, ctx.get_stream<cpu>());
+ } else {
+ std::vector<AType> a_vec(a.Size(), 0);
+ std::vector<WType> w_vec(w.Size(), 0);
+ AType* a_cp_ptr = a_vec.data();
+ WType* w_cp_ptr = w_vec.data();
+ GpuCallbackCpuImpl(a, w, a_cp_ptr, w_cp_ptr, &workspace, ctx, req);
+ }
+ });
+ });
+}
+
+struct EigvalshParam : public dmlc::Parameter<EigvalshParam> {
+ char UPLO;
+ DMLC_DECLARE_PARAMETER(EigvalshParam) {
+ DMLC_DECLARE_FIELD(UPLO)
+ .set_default('L')
+ .describe("Specifies whether the calculation is done with the lower or upper triangular part.");
+ }
+};
+
+template<typename xpu>
+size_t EighEigvalshForwardWorkspaceSize(const TBlob& a,
+ const TBlob& w,
Review comment:
alignment
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364998303
##########
File path: src/operator/numpy/linalg/np_eig.cu
##########
@@ -0,0 +1,61 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * \file np_eigvals.cu
+ * \brief GPU implementation placeholder of Eigvals Operator
+ */
+
+#include <mxnet/operator_util.h>
+#include "./np_eig-inl.h"
+
+namespace mxnet {
+namespace op {
+
+template<>
+size_t EigForwardWorkspaceSize<gpu>(const TBlob& a,
+ const TBlob& w,
+ const TBlob& v,
+ const std::vector<OpReqType>& req) {
+ LOG(FATAL) << "EigForwardWorkspaceSize in gpu is unsupported";
+ return 0U;
+}
+
+template<>
+void EigOpForwardImpl<gpu>(const TBlob& a,
Review comment:
Those 2 functions are never called, you could probably delete them.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364996931
##########
File path: src/operator/numpy/linalg/np_eigvals.cc
##########
@@ -0,0 +1,124 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals.cc
+ * \brief CPU implementation placeholder of Eigvals Operator
+ */
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig.
+bool EigvalsOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape)) {
Review comment:
you could probably use
```c++
if (shape_is_known(a_shape)) {
// code
} else if (shape_is_known(eig_shape)) {
// code
}
```
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r366024841
##########
File path: src/operator/numpy/linalg/np_eigvals-inl.h
##########
@@ -0,0 +1,432 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals-inl.h
+ * \brief Placeholder for eigvals
+ */
+#ifndef MXNET_OPERATOR_NUMPY_LINALG_NP_EIGVALS_INL_H_
+#define MXNET_OPERATOR_NUMPY_LINALG_NP_EIGVALS_INL_H_
+
+#include <mxnet/operator_util.h>
+#include <vector>
+#include "../../operator_common.h"
+#include "../../mshadow_op.h"
+#include "../../tensor/la_op.h"
+#include "../../tensor/la_op-inl.h"
+#include "./np_solve-inl.h"
+
+namespace mxnet {
+namespace op {
+
+using namespace mshadow;
+
+template<int req>
+struct eigvals_assign_helper {
+ template<typename DType>
+ MSHADOW_XINLINE static void Map(int i, const DType *in_data, DType *out_data) {
+ KERNEL_ASSIGN(out_data[i], req, in_data[i]);
+ }
+};
+
+struct eigh_eigvalsh_helper {
+ template<typename InDType, typename OutDType>
+ MSHADOW_XINLINE static void Map(int i, const InDType *in_data, OutDType *out_data,
+ const int nrow, const int ld, const int step, bool USE_UP) {
+ int idx = i / step, row = (i % step) / ld, col = (i % step) % ld;
+ if ((USE_UP && row > col) || (!USE_UP && row < col)) {
+ out_data[idx * step + col + row * ld] =
+ static_cast<OutDType>(in_data[idx * step + row + col * ld]);
+ } else {
+ out_data[idx * step + col + row * ld] =
+ static_cast<OutDType>(in_data[idx * step + col + row * ld]);
+ }
+ }
+};
+
+template<typename xpu, typename DType>
+void linalg_geev(char jobvl,
+ char jobvr,
+ const Tensor<xpu, 2, DType>& a,
+ const Tensor<xpu, 1, DType>& wr,
+ const Tensor<xpu, 1, DType>& wi,
+ const Tensor<xpu, 2, DType>& vl,
+ const Tensor<xpu, 2, DType>& vr,
+ const Tensor<xpu, 1, DType>& work_array);
+
+#define LINALG_CPU_EIG(fname, DType) \
+template<> inline \
+void linalg_geev<cpu, DType>(char jobvl, \
+ char jobvr, \
+ const Tensor<cpu, 2, DType>& a, \
+ const Tensor<cpu, 1, DType>& wr, \
+ const Tensor<cpu, 1, DType>& wi, \
+ const Tensor<cpu, 2, DType>& vl, \
+ const Tensor<cpu, 2, DType>& vr, \
+ const Tensor<cpu, 1, DType>& work_array) { \
+ const int n = a.size(1), lda = a.size(0); \
+ const int lwork = work_array.shape_.Size(); \
+ const int ldvl = vl.size(0), ldvr = vr.size(0); \
+ int res(MXNET_LAPACK_##fname(MXNET_LAPACK_COL_MAJOR, jobvl, jobvr, \
+ n, a.dptr_, lda, \
+ wr.dptr_, wi.dptr_, \
+ vl.dptr_, ldvl, \
+ vr.dptr_, ldvr, \
+ work_array.dptr_, lwork)); \
+ CHECK_LE(res, 0) << #fname << "the QR algorithm failed to compute all the" \
+ << "eigenvalues, and no eigenvectors have been computed; elements " \
+ << res + 1 << ":N" << " of WR and WI contain eigenvalues which have converged"; \
+ CHECK_GE(res, 0) << #fname << ": the " << -res \
+ << "-th argument had an illegal value"; \
+}
Review comment:
blank line below.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r364461515
##########
File path: src/operator/numpy/linalg/np_eig-inl.h
##########
@@ -0,0 +1,274 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eig-inl.h
+ * \brief Placeholder for eig
+ */
+#ifndef MXNET_OPERATOR_NUMPY_LINALG_NP_EIG_INL_H_
+#define MXNET_OPERATOR_NUMPY_LINALG_NP_EIG_INL_H_
+
+#include <vector>
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+using namespace mshadow;
+
+template<int req>
+struct eigvec_assign_helper {
+ template<typename DType>
+ MSHADOW_XINLINE static void Map(int i, const DType *in_data, DType *out_data,
+ const int nrow, const int ld, const int step) {
+ int idx = i / step, row = (i % step) / ld, col = (i % step) % ld;
+ KERNEL_ASSIGN(out_data[idx * step + row + col * ld], req, in_data[i]);
+ }
+};
+
+// Calculates workspace size of eig forward op.
+// The dimension of the array WORK in LAPACKE_#GEEV should >= max(1,3*N), and
+// if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.
+// For good performance, LWORK must generally be larger.
+template<typename xpu>
+size_t EigForwardWorkspaceSize(const TBlob& a,
+ const TBlob& w,
+ const TBlob& v,
+ const std::vector<OpReqType>& req) {
+ if (kNullOp == req[0] && kNullOp == req[1]) { return 0U; }
+
+ // Zero-size input, no need to launch kernel
+ if (0U == a.Size()) { return 0U; }
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, DType, {
+ size_t work_space_size = 0;
+ size_t n = a.size(a.ndim() - 1);
+ work_space_size += a.Size(); // For matrix.
+ work_space_size += 2 * w.Size(); // For eigenvalues' real and image component.
+ work_space_size += n * n; // For left eigenvectors temp memory
+ work_space_size += v.Size(); // For right eigenvectors real and image component.
+ work_space_size += 4 * n; // For workspace size in LAPACKE_#GEEV.
+ work_space_size *= sizeof(DType);
+ return work_space_size;
+ });
+ LOG(FATAL) << "InternalError: cannot reach here";
+ return 0U;
+}
+
+template<typename xpu>
+void EigOpForwardImpl(const TBlob& a,
+ const TBlob& w,
+ const TBlob& v,
+ const std::vector<OpReqType>& req,
+ std::vector<char> *workspace,
+ mshadow::Stream<xpu> *s) {
+ if (kNullOp == req[0] && kNullOp == req[1]) { return; }
+ const mxnet::TShape& a_shape = a.shape_;
+ const mxnet::TShape& w_shape = w.shape_;
+ const mxnet::TShape& v_shape = v.shape_;
+ const int a_ndim = a_shape.ndim();
+
+ // Zero-size output, no need to launch kernel
+ if (0U == a.Size()) { return; }
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, DType, {
+ const int N = a_shape[a_ndim - 1];
+ DType *a_ptr =
+ reinterpret_cast<DType*>(workspace->data());
+ DType *wr_ptr =
+ reinterpret_cast<DType*>(workspace->data() + a.Size() * sizeof(DType));
+ DType *wi_ptr =
+ reinterpret_cast<DType*>(workspace->data() + (w.Size() + a.Size()) * sizeof(DType));
+ DType *vl_ptr =
+ reinterpret_cast<DType*>(workspace->data() + (2 * w.Size() + a.Size()) * sizeof(DType));
+ DType *vr_ptr =
+ reinterpret_cast<DType*>(
+ workspace->data() + (2 * w.Size() + N * N + a.Size()) * sizeof(DType));
+ DType *work_ptr =
+ reinterpret_cast<DType*>(
+ workspace->data() + (2 * w.Size() + v.Size() + N * N + a.Size()) * sizeof(DType));
+ MSHADOW_TYPE_SWITCH(a.type_flag_, AType, {
+ // Cast type and transpose.
+ mxnet_op::Kernel<SolveTypeTransposeHelper, xpu>::Launch(
+ s, a_shape.Size(), a.dptr<AType>(), a_ptr, N, N, N * N);
+ });
+ char jobvl = 'N', jobvr = 'V';
+ mxnet::TBlob a_trans_data(a_ptr, a_shape, a.dev_mask(), a.dev_id());
+ mxnet::TBlob wr_data(wr_ptr, w_shape, w.dev_mask(), w.dev_id());
+ mxnet::TBlob wi_data(wi_ptr, w_shape, w.dev_mask(), w.dev_id());
+ mxnet::TBlob vl_data(vl_ptr, Shape3(1, N, N), v.dev_mask(), v.dev_id());
+ mxnet::TBlob vr_data(vr_ptr, v_shape, v.dev_mask(), v.dev_id());
+ mxnet::TBlob work_data(work_ptr, Shape1(4 * N), a.dev_mask(), a.dev_id());
+ eig_eigvals::op(jobvl, jobvr,
+ a_trans_data.FlatToKD<xpu, 3, DType>(s),
+ wr_data.FlatToKD<xpu, 2, DType>(s),
+ wi_data.FlatToKD<xpu, 2, DType>(s),
+ vl_data.get<xpu, 3, DType>(s),
+ vr_data.FlatToKD<xpu, 3, DType>(s),
+ work_data.get<xpu, 1, DType>(s));
+ for (size_t i = 0; i < wi_data.Size(); ++i) {
+ CHECK_LE(fabs(wi_ptr[i]), 1e-15)
+ << "Complex eigvals is unsupported in linalg temporary.";
+ }
+ MXNET_ASSIGN_REQ_SWITCH(req[0], req_type, {
+ mxnet_op::Kernel<eigvals_assign_helper<req_type>, xpu>::Launch(
+ s, w.Size(), wr_ptr, w.dptr<DType>());
+ });
+ MXNET_ASSIGN_REQ_SWITCH(req[1], req_type, {
+ mxnet_op::Kernel<eigvec_assign_helper<req_type>, xpu>::Launch(
+ s, v.Size(), vr_ptr, v.dptr<DType>(), N, N, N * N);
+ });
+ });
+}
+
+template<typename AType, typename WType>
+void GpuCallbackCpuImpl(const TBlob& a,
+ const TBlob& w,
+ const TBlob& v,
+ AType* a_cp_ptr,
+ WType* w_cp_ptr,
+ WType* v_cp_ptr,
+ std::vector<char>* workspace,
+ const OpContext& ctx,
+ const std::vector<OpReqType>& req) {
+#if MXNET_USE_CUDA
+ mshadow::Stream<gpu> *s = ctx.get_stream<gpu>();
+ cudaStream_t stream = Stream<gpu>::GetStream(s);
+ CUDA_CALL(cudaMemcpyAsync(a_cp_ptr, a.dptr<AType>(), sizeof(AType) * a.Size(),
+ cudaMemcpyDeviceToHost, stream));
+ CUDA_CALL(cudaMemcpyAsync(w_cp_ptr, w.dptr<WType>(), sizeof(WType) * w.Size(),
+ cudaMemcpyDeviceToHost, stream));
+ CUDA_CALL(cudaMemcpyAsync(v_cp_ptr, v.dptr<WType>(), sizeof(WType) * v.Size(),
+ cudaMemcpyDeviceToHost, stream));
+ CUDA_CALL(cudaStreamSynchronize(stream));
+ mxnet::TBlob a_data(a_cp_ptr, a.shape_, cpu::kDevMask);
+ mxnet::TBlob w_data(w_cp_ptr, w.shape_, cpu::kDevMask);
+ mxnet::TBlob v_data(v_cp_ptr, v.shape_, cpu::kDevMask);
+ // Op forward implement on cpu.
+ EigOpForwardImpl<cpu>(a_data, w_data, v_data, req, workspace, ctx.get_stream<cpu>());
+ // Copy back to gpu.
+ CUDA_CALL(cudaMemcpyAsync(w.dptr<WType>(), w_cp_ptr, sizeof(WType) * w.Size(),
+ cudaMemcpyHostToDevice, stream));
+ CUDA_CALL(cudaMemcpyAsync(v.dptr<WType>(), v_cp_ptr, sizeof(WType) * v.Size(),
+ cudaMemcpyHostToDevice, stream));
+ CUDA_CALL(cudaStreamSynchronize(stream));
+#else
+ LOG(FATAL) << "Please build with USE_CUDA=1 to enable GPU";
+#endif // MXNET_USE_CUDA
+}
+
+template<typename xpu>
+void EigOpForward(const nnvm::NodeAttrs& attrs,
+ const OpContext& ctx,
+ const std::vector<TBlob>& inputs,
+ const std::vector<OpReqType>& req,
+ const std::vector<TBlob>& outputs) {
+ CHECK_EQ(inputs.size(), 1U);
+ CHECK_EQ(outputs.size(), 2U);
+ CHECK_EQ(req.size(), 2U);
+ const TBlob& a = inputs[0];
+ const TBlob& w = outputs[0];
+ const TBlob& v = outputs[1];
+
+ // Calculate workspace size.
+ size_t workspace_size = EigForwardWorkspaceSize<cpu>(a, w, v, req);
+ std::vector<char> workspace(workspace_size, 0);
+
+ MSHADOW_SGL_DBL_TYPE_SWITCH(w.type_flag_, WType, {
+ MSHADOW_SGL_DBL_TYPE_SWITCH(a.type_flag_, AType, {
+ if (xpu::kDevCPU) {
+ // Op forward implement.
+ EigOpForwardImpl<cpu>(a, w, v, req, &workspace, ctx.get_stream<cpu>());
+ } else {
+ std::vector<AType> a_vec(a.Size(), 0);
+ std::vector<WType> w_vec(w.Size(), 0);
+ std::vector<WType> v_vec(v.Size(), 0);
+ AType* a_cp_ptr = a_vec.data();
+ WType* w_cp_ptr = w_vec.data();
+ WType* v_cp_ptr = v_vec.data();
+ GpuCallbackCpuImpl(a, w, v, a_cp_ptr, w_cp_ptr, v_cp_ptr, &workspace, ctx, req);
Review comment:
suggested:
```c++
#if MXNET_USE_CUDA
template<typename AType, typename WType>
void GpuCallbackCpuImpl(...) {
}
#endif
// ...
if (xpu::kDevCPU) {
// ...
} else {
#if MXNET_USE_CUDA
// ...
GpuCallbackCpuImpl(...);
#else
LOG(FATAL) << "Please build with USE_CUDA=1 to enable GPU";
#endif
}
```
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r365002633
##########
File path: src/operator/numpy/linalg/np_eig.cc
##########
@@ -0,0 +1,157 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eig.cc
+ * \brief CPU implementation placeholder of Eig Operator
+ */
+#include "./np_eig-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig, EigVec
+bool EigOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 2U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+ const mxnet::TShape& eigv_shape = (*out_attrs)[1];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ // Calculate eig vec shape: must have the same shape as A
+ SHAPE_ASSIGN_CHECK(*out_attrs, 1, a_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape) && shape_is_known(eigv_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ const int eigv_ndim = eigv_shape.ndim();
+ CHECK_GE(eigv_ndim, 2)
+ << "Outputs V must be at least two-dimensional";
+ CHECK_EQ(eigv_shape[eigv_ndim - 2], eigv_shape[eigv_ndim - 1])
+ << "Outputs V's last two dimension must be equal";
+ CHECK_EQ(eig_ndim + 1, eigv_ndim)
+ << "Outputs W, V must satisfy W.ndim == V.ndim - 1";
+ for (int i = 0; i < eig_ndim; ++i) {
+ CHECK_EQ(eig_shape[i], eigv_shape[i])
+ << "Outputs W, V's shape dismatch";
+ }
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, eigv_shape);
+ } else if (shape_is_known(eig_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ CHECK_GE(eig_ndim, 1)
+ << "Outputs W must be at least one-dimensional";
+ std::vector<int> eigv_shape_vec(eig_ndim + 1);
+ for (int i = 0; i < eig_ndim; ++i) {
+ eigv_shape_vec[i] = eig_shape[i];
+ }
+ eigv_shape_vec[eig_ndim] = eig_shape[eig_ndim - 1];
+ mxnet::TShape eigv_shape(eigv_shape_vec.begin(), eigv_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, eigv_shape);
+ SHAPE_ASSIGN_CHECK(*out_attrs, 1, eigv_shape);
+ } else {
+ const int eigv_ndim = eigv_shape.ndim();
+ CHECK_GE(eigv_ndim, 2)
+ << "Outputs V must be at least two-dimensional";
+ CHECK_EQ(eigv_shape[eigv_ndim - 2], eigv_shape[eigv_ndim - 1])
+ << "Outputs V's last two dimension must be equal";
+ std::vector<int> eig_shape_vec(eigv_ndim - 1);
+ for (int i = 0; i < eigv_ndim - 1; ++i) {
+ eig_shape_vec[i] = eigv_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, eigv_shape);
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ }
+ }
+ return shape_is_known(*in_attrs) && shape_is_known(*out_attrs);
+}
+
+inline bool EigOpType(const nnvm::NodeAttrs& attrs,
+ std::vector<int>* in_attrs,
+ std::vector<int>* out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
Review comment:
see the comment below on the other InferType function.
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services
[GitHub] [incubator-mxnet] haojin2 commented on a change in pull request
#17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
Posted by GitBox <gi...@apache.org>.
haojin2 commented on a change in pull request #17188: [Numpy] Add linalg.eig/eigh/eigvals/eigvalsh op
URL: https://github.com/apache/incubator-mxnet/pull/17188#discussion_r365001641
##########
File path: src/operator/numpy/linalg/np_eigvals.cc
##########
@@ -0,0 +1,124 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://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.
+ */
+
+/*!
+ * Copyright (c) 2019 by Contributors
+ * \file np_eigvals.cc
+ * \brief CPU implementation placeholder of Eigvals Operator
+ */
+#include "./np_eigvals-inl.h"
+
+namespace mxnet {
+namespace op {
+
+// Inputs: A.
+// Outputs: Eig.
+bool EigvalsOpShape(const nnvm::NodeAttrs& attrs,
+ mxnet::ShapeVector *in_attrs,
+ mxnet::ShapeVector *out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ const mxnet::TShape& a_shape = (*in_attrs)[0];
+ const mxnet::TShape& eig_shape = (*out_attrs)[0];
+
+ if (shape_is_known(a_shape)) {
+ // Forward shape inference.
+ const int a_ndim = a_shape.ndim();
+ CHECK_GE(a_ndim, 2)
+ << "Array must be at least two-dimensional";
+ CHECK_EQ(a_shape[a_ndim - 2], a_shape[a_ndim - 1])
+ << "Input A's last two dimension must be equal";
+
+ // Calculate eig shape.
+ std::vector<int> eig_shape_vec(a_ndim - 1, -1);
+ for (int i = 0; i < a_ndim - 1; ++i) {
+ eig_shape_vec[i] = a_shape[i];
+ }
+ mxnet::TShape eig_shape(eig_shape_vec.begin(), eig_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*out_attrs, 0, eig_shape);
+ } else {
+ // Backward shape inference.
+ if (shape_is_known(eig_shape)) {
+ const int eig_ndim = eig_shape.ndim();
+ CHECK_GE(eig_ndim, 1)
+ << "Outputs W must be at least one-dimensional";
+ std::vector<int> a_shape_vec(eig_ndim + 1);
+ for (int i = 0; i < eig_ndim; ++i) {
+ a_shape_vec[i] = eig_shape[i];
+ }
+ a_shape_vec[eig_ndim] = eig_shape[eig_ndim - 1];
+ mxnet::TShape a_shape(a_shape_vec.begin(), a_shape_vec.end());
+ SHAPE_ASSIGN_CHECK(*in_attrs, 0, a_shape);
+ }
+ }
+ return shape_is_known(*in_attrs) && shape_is_known(*out_attrs);
+}
+
+inline bool EigvalsOpType(const nnvm::NodeAttrs& attrs,
+ std::vector<int>* in_attrs,
+ std::vector<int>* out_attrs) {
+ CHECK_EQ(in_attrs->size(), 1U);
+ CHECK_EQ(out_attrs->size(), 1U);
+ int a_type = in_attrs->at(0);
+ // unsupport float16
Review comment:
Does the 2 ops support integer inputs?
If not, then I would suggest the following code structure:
```c++
if (common::is_float(a_type)) {
// raise error for fp16
TYPE_ASSIGN_CHECK(*out_attrs, 0, in_attrs->at(0));
} else {
// raise error
}
```
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
users@infra.apache.org
With regards,
Apache Git Services