We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
`cupy.linalg.solve` returns wrong output for complex inputs
```
>>> a
array([[ 1.03310883-0.87701077j, -0.35775642-0.94918977j,
         2.63837642+0.29920745j],
       [-0.88041189+0.97322198j,  0.26771522-0.56788898j,
        -0.47020123+0.63015159j],
       [-0.05746688-0.16238003j, -0.18490837-0.00778188j,
         1.61550173+1.13397892j]])
>>> b
array([ 0.92632861+1.45552766j, -0.54307601-0.62108012j,
       -0.23769241+2.91680425j])
>>> cupy.linalg.solve(a, b)
array([ 0.07332535+0.27916782j, -1.61809337-0.62108012j,
       -0.23769241+2.91680425j])
>>> numpy.linalg.solve(cupy.asnumpy(a), cupy.asnumpy(b))
array([0.04332182-0.34957548j, 1.6260605 -0.2383031j ,
       0.88753974+1.15498982j])
>>> cupy.__version__
'5.0.0b4'
```

#1524 would fix the issue.

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `cupy/linalg/solve.py`
Content:
```
1 import numpy
2 from numpy import linalg
3 import six
4 
5 import cupy
6 from cupy.core import core
7 from cupy import cuda
8 from cupy.cuda import cublas
9 from cupy.cuda import device
10 from cupy.linalg import decomposition
11 from cupy.linalg import util
12 
13 if cuda.cusolver_enabled:
14     from cupy.cuda import cusolver
15 
16 
17 def solve(a, b):
18     """Solves a linear matrix equation.
19 
20     It computes the exact solution of ``x`` in ``ax = b``,
21     where ``a`` is a square and full rank matrix.
22 
23     Args:
24         a (cupy.ndarray): The matrix with dimension ``(..., M, M)``.
25         b (cupy.ndarray): The matrix with dimension ``(...,M)`` or
26             ``(..., M, K)``.
27 
28     Returns:
29         cupy.ndarray:
30             The matrix with dimension ``(..., M)`` or ``(..., M, K)``.
31 
32     .. seealso:: :func:`numpy.linalg.solve`
33     """
34     # NOTE: Since cusolver in CUDA 8.0 does not support gesv,
35     #       we manually solve a linear system with QR decomposition.
36     #       For details, please see the following:
37     #       https://docs.nvidia.com/cuda/cusolver/index.html#qr_examples
38     if not cuda.cusolver_enabled:
39         raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')
40 
41     util._assert_cupy_array(a, b)
42     util._assert_nd_squareness(a)
43 
44     if not ((a.ndim == b.ndim or a.ndim == b.ndim + 1) and
45             a.shape[:-1] == b.shape[:a.ndim - 1]):
46         raise ValueError(
47             'a must have (..., M, M) shape and b must have (..., M) '
48             'or (..., M, K)')
49 
50     # Cast to float32 or float64
51     if a.dtype.char == 'f' or a.dtype.char == 'd':
52         dtype = a.dtype
53     else:
54         dtype = numpy.find_common_type((a.dtype.char, 'f'), ())
55 
56     a = a.astype(dtype)
57     b = b.astype(dtype)
58     if a.ndim == 2:
59         return _solve(a, b)
60     x = cupy.empty_like(b)
61     shape = a.shape[:-2]
62     for i in six.moves.range(numpy.prod(shape)):
63         index = numpy.unravel_index(i, shape)
64         x[index] = _solve(a[index], b[index])
65     return x
66 
67 
68 def _solve(a, b):
69     a = cupy.asfortranarray(a)
70     b = cupy.asfortranarray(b)
71     dtype = a.dtype
72     m, k = (b.size, 1) if b.ndim == 1 else b.shape
73     cusolver_handle = device.get_cusolver_handle()
74     cublas_handle = device.get_cublas_handle()
75     dev_info = cupy.empty(1, dtype=numpy.int32)
76 
77     if dtype == 'f':
78         geqrf = cusolver.sgeqrf
79         geqrf_bufferSize = cusolver.sgeqrf_bufferSize
80         ormqr = cusolver.sormqr
81         trsm = cublas.strsm
82     else:  # dtype == 'd'
83         geqrf = cusolver.dgeqrf
84         geqrf_bufferSize = cusolver.dgeqrf_bufferSize
85         ormqr = cusolver.dormqr
86         trsm = cublas.dtrsm
87 
88     # 1. QR decomposition (A = Q * R)
89     buffersize = geqrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)
90     workspace = cupy.empty(buffersize, dtype=dtype)
91     tau = cupy.empty(m, dtype=dtype)
92     geqrf(
93         cusolver_handle, m, m, a.data.ptr, m,
94         tau.data.ptr, workspace.data.ptr, buffersize, dev_info.data.ptr)
95     _check_status(dev_info)
96     # 2. ormqr (Q^T * B)
97     ormqr(
98         cusolver_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_OP_T,
99         m, k, m, a.data.ptr, m, tau.data.ptr, b.data.ptr, m,
100         workspace.data.ptr, buffersize, dev_info.data.ptr)
101     _check_status(dev_info)
102     # 3. trsm (X = R^{-1} * (Q^T * B))
103     trsm(
104         cublas_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_FILL_MODE_UPPER,
105         cublas.CUBLAS_OP_N, cublas.CUBLAS_DIAG_NON_UNIT,
106         m, k, 1, a.data.ptr, m, b.data.ptr, m)
107     return b
108 
109 
110 def _check_status(dev_info):
111     status = int(dev_info)
112     if status < 0:
113         raise linalg.LinAlgError(
114             'Parameter error (maybe caused by a bug in cupy.linalg?)')
115 
116 
117 def tensorsolve(a, b, axes=None):
118     """Solves tensor equations denoted by ``ax = b``.
119 
120     Suppose that ``b`` is equivalent to ``cupy.tensordot(a, x)``.
121     This function computes tensor ``x`` from ``a`` and ``b``.
122 
123     Args:
124         a (cupy.ndarray): The tensor with ``len(shape) >= 1``
125         b (cupy.ndarray): The tensor with ``len(shape) >= 1``
126         axes (tuple of ints): Axes in ``a`` to reorder to the right
127             before inversion.
128 
129     Returns:
130         cupy.ndarray:
131             The tensor with shape ``Q`` such that ``b.shape + Q == a.shape``.
132 
133     .. seealso:: :func:`numpy.linalg.tensorsolve`
134     """
135     if axes is not None:
136         allaxes = list(six.moves.range(a.ndim))
137         for k in axes:
138             allaxes.remove(k)
139             allaxes.insert(a.ndim, k)
140         a = a.transpose(allaxes)
141 
142     oldshape = a.shape[-(a.ndim - b.ndim):]
143     prod = cupy.internal.prod(oldshape)
144 
145     a = a.reshape(-1, prod)
146     b = b.ravel()
147     result = solve(a, b)
148     return result.reshape(oldshape)
149 
150 
151 # TODO(okuta): Implement lstsq
152 
153 
154 def inv(a):
155     """Computes the inverse of a matrix.
156 
157     This function computes matrix ``a_inv`` from n-dimensional regular matrix
158     ``a`` such that ``dot(a, a_inv) == eye(n)``.
159 
160     Args:
161         a (cupy.ndarray): The regular matrix
162 
163     Returns:
164         cupy.ndarray: The inverse of a matrix.
165 
166     .. seealso:: :func:`numpy.linalg.inv`
167     """
168     if not cuda.cusolver_enabled:
169         raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')
170 
171     # to prevent `a` to be overwritten
172     a = a.copy()
173 
174     util._assert_cupy_array(a)
175     util._assert_rank2(a)
176     util._assert_nd_squareness(a)
177 
178     if a.dtype.char == 'f' or a.dtype.char == 'd':
179         dtype = a.dtype.char
180     else:
181         dtype = numpy.find_common_type((a.dtype.char, 'f'), ()).char
182 
183     cusolver_handle = device.get_cusolver_handle()
184     dev_info = cupy.empty(1, dtype=dtype)
185 
186     ipiv = cupy.empty((a.shape[0], 1), dtype=dtype)
187 
188     if dtype == 'f':
189         getrf = cusolver.sgetrf
190         getrf_bufferSize = cusolver.sgetrf_bufferSize
191         getrs = cusolver.sgetrs
192     else:  # dtype == 'd'
193         getrf = cusolver.dgetrf
194         getrf_bufferSize = cusolver.dgetrf_bufferSize
195         getrs = cusolver.dgetrs
196 
197     m = a.shape[0]
198 
199     buffersize = getrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)
200     workspace = cupy.empty(buffersize, dtype=dtype)
201 
202     # LU factorization
203     getrf(cusolver_handle, m, m, a.data.ptr, m, workspace.data.ptr,
204           ipiv.data.ptr, dev_info.data.ptr)
205 
206     b = cupy.eye(m, dtype=dtype)
207 
208     # solve for the inverse
209     getrs(cusolver_handle, 0, m, m, a.data.ptr, m, ipiv.data.ptr, b.data.ptr,
210           m, dev_info.data.ptr)
211 
212     return b
213 
214 
215 def pinv(a, rcond=1e-15):
216     """Compute the Moore-Penrose pseudoinverse of a matrix.
217 
218     It computes a pseudoinverse of a matrix ``a``, which is a generalization
219     of the inverse matrix with Singular Value Decomposition (SVD).
220     Note that it automatically removes small singular values for stability.
221 
222     Args:
223         a (cupy.ndarray): The matrix with dimension ``(M, N)``
224         rcond (float): Cutoff parameter for small singular values.
225             For stability it computes the largest singular value denoted by
226             ``s``, and sets all singular values smaller than ``s`` to zero.
227 
228     Returns:
229         cupy.ndarray: The pseudoinverse of ``a`` with dimension ``(N, M)``.
230 
231     .. seealso:: :func:`numpy.linalg.pinv`
232     """
233     u, s, vt = decomposition.svd(a, full_matrices=False)
234     cutoff = rcond * s.max()
235     s1 = 1 / s
236     s1[s <= cutoff] = 0
237     return core.dot(vt.T, s1[:, None] * u.T)
238 
239 
240 def tensorinv(a, ind=2):
241     """Computes the inverse of a tensor.
242 
243     This function computes tensor ``a_inv`` from tensor ``a`` such that
244     ``tensordot(a_inv, a, ind) == I``, where ``I`` denotes the identity tensor.
245 
246     Args:
247         a (cupy.ndarray):
248             The tensor such that
249             ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
250         ind (int):
251             The positive number used in ``axes`` option of ``tensordot``.
252 
253     Returns:
254         cupy.ndarray:
255             The inverse of a tensor whose shape is equivalent to
256             ``a.shape[ind:] + a.shape[:ind]``.
257 
258     .. seealso:: :func:`numpy.linalg.tensorinv`
259     """
260     util._assert_cupy_array(a)
261 
262     if ind <= 0:
263         raise ValueError('Invalid ind argument')
264     oldshape = a.shape
265     invshape = oldshape[ind:] + oldshape[:ind]
266     prod = cupy.internal.prod(oldshape[ind:])
267     a = a.reshape(prod, -1)
268     a_inv = inv(a)
269     return a_inv.reshape(*invshape)
270
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/cupy/linalg/solve.py b/cupy/linalg/solve.py
--- a/cupy/linalg/solve.py
+++ b/cupy/linalg/solve.py
@@ -53,37 +53,55 @@
     else:
         dtype = numpy.find_common_type((a.dtype.char, 'f'), ())
 
+    cublas_handle = device.get_cublas_handle()
+    cusolver_handle = device.get_cusolver_handle()
+
     a = a.astype(dtype)
     b = b.astype(dtype)
     if a.ndim == 2:
-        return _solve(a, b)
+        return _solve(a, b, cublas_handle, cusolver_handle)
+
     x = cupy.empty_like(b)
     shape = a.shape[:-2]
     for i in six.moves.range(numpy.prod(shape)):
         index = numpy.unravel_index(i, shape)
-        x[index] = _solve(a[index], b[index])
+        x[index] = _solve(a[index], b[index], cublas_handle, cusolver_handle)
     return x
 
 
-def _solve(a, b):
+def _solve(a, b, cublas_handle, cusolver_handle):
     a = cupy.asfortranarray(a)
     b = cupy.asfortranarray(b)
     dtype = a.dtype
     m, k = (b.size, 1) if b.ndim == 1 else b.shape
-    cusolver_handle = device.get_cusolver_handle()
-    cublas_handle = device.get_cublas_handle()
     dev_info = cupy.empty(1, dtype=numpy.int32)
 
     if dtype == 'f':
         geqrf = cusolver.sgeqrf
         geqrf_bufferSize = cusolver.sgeqrf_bufferSize
         ormqr = cusolver.sormqr
+        trans = cublas.CUBLAS_OP_T
         trsm = cublas.strsm
-    else:  # dtype == 'd'
+    elif dtype == 'd':
         geqrf = cusolver.dgeqrf
         geqrf_bufferSize = cusolver.dgeqrf_bufferSize
         ormqr = cusolver.dormqr
+        trans = cublas.CUBLAS_OP_T
         trsm = cublas.dtrsm
+    elif dtype == 'F':
+        geqrf = cusolver.cgeqrf
+        geqrf_bufferSize = cusolver.cgeqrf_bufferSize
+        ormqr = cusolver.cormqr
+        trans = cublas.CUBLAS_OP_C
+        trsm = cublas.ctrsm
+    elif dtype == 'D':
+        geqrf = cusolver.zgeqrf
+        geqrf_bufferSize = cusolver.zgeqrf_bufferSize
+        ormqr = cusolver.zormqr
+        trans = cublas.CUBLAS_OP_C
+        trsm = cublas.ztrsm
+    else:
+        raise NotImplementedError(dtype)
 
     # 1. QR decomposition (A = Q * R)
     buffersize = geqrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)
@@ -95,7 +113,7 @@
     _check_status(dev_info)
     # 2. ormqr (Q^T * B)
     ormqr(
-        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_OP_T,
+        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, trans,
         m, k, m, a.data.ptr, m, tau.data.ptr, b.data.ptr, m,
         workspace.data.ptr, buffersize, dev_info.data.ptr)
     _check_status(dev_info)

{"golden_diff": "diff --git a/cupy/linalg/solve.py b/cupy/linalg/solve.py\n--- a/cupy/linalg/solve.py\n+++ b/cupy/linalg/solve.py\n@@ -53,37 +53,55 @@\n     else:\n         dtype = numpy.find_common_type((a.dtype.char, 'f'), ())\n \n+    cublas_handle = device.get_cublas_handle()\n+    cusolver_handle = device.get_cusolver_handle()\n+\n     a = a.astype(dtype)\n     b = b.astype(dtype)\n     if a.ndim == 2:\n-        return _solve(a, b)\n+        return _solve(a, b, cublas_handle, cusolver_handle)\n+\n     x = cupy.empty_like(b)\n     shape = a.shape[:-2]\n     for i in six.moves.range(numpy.prod(shape)):\n         index = numpy.unravel_index(i, shape)\n-        x[index] = _solve(a[index], b[index])\n+        x[index] = _solve(a[index], b[index], cublas_handle, cusolver_handle)\n     return x\n \n \n-def _solve(a, b):\n+def _solve(a, b, cublas_handle, cusolver_handle):\n     a = cupy.asfortranarray(a)\n     b = cupy.asfortranarray(b)\n     dtype = a.dtype\n     m, k = (b.size, 1) if b.ndim == 1 else b.shape\n-    cusolver_handle = device.get_cusolver_handle()\n-    cublas_handle = device.get_cublas_handle()\n     dev_info = cupy.empty(1, dtype=numpy.int32)\n \n     if dtype == 'f':\n         geqrf = cusolver.sgeqrf\n         geqrf_bufferSize = cusolver.sgeqrf_bufferSize\n         ormqr = cusolver.sormqr\n+        trans = cublas.CUBLAS_OP_T\n         trsm = cublas.strsm\n-    else:  # dtype == 'd'\n+    elif dtype == 'd':\n         geqrf = cusolver.dgeqrf\n         geqrf_bufferSize = cusolver.dgeqrf_bufferSize\n         ormqr = cusolver.dormqr\n+        trans = cublas.CUBLAS_OP_T\n         trsm = cublas.dtrsm\n+    elif dtype == 'F':\n+        geqrf = cusolver.cgeqrf\n+        geqrf_bufferSize = cusolver.cgeqrf_bufferSize\n+        ormqr = cusolver.cormqr\n+        trans = cublas.CUBLAS_OP_C\n+        trsm = cublas.ctrsm\n+    elif dtype == 'D':\n+        geqrf = cusolver.zgeqrf\n+        geqrf_bufferSize = cusolver.zgeqrf_bufferSize\n+        ormqr = cusolver.zormqr\n+        trans = cublas.CUBLAS_OP_C\n+        trsm = cublas.ztrsm\n+    else:\n+        raise NotImplementedError(dtype)\n \n     # 1. QR decomposition (A = Q * R)\n     buffersize = geqrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)\n@@ -95,7 +113,7 @@\n     _check_status(dev_info)\n     # 2. ormqr (Q^T * B)\n     ormqr(\n-        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_OP_T,\n+        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, trans,\n         m, k, m, a.data.ptr, m, tau.data.ptr, b.data.ptr, m,\n         workspace.data.ptr, buffersize, dev_info.data.ptr)\n     _check_status(dev_info)\n", "issue": "`cupy.linalg.solve` returns wrong output for complex inputs\n```\r\n>>> a\r\narray([[ 1.03310883-0.87701077j, -0.35775642-0.94918977j,\r\n         2.63837642+0.29920745j],\r\n       [-0.88041189+0.97322198j,  0.26771522-0.56788898j,\r\n        -0.47020123+0.63015159j],\r\n       [-0.05746688-0.16238003j, -0.18490837-0.00778188j,\r\n         1.61550173+1.13397892j]])\r\n>>> b\r\narray([ 0.92632861+1.45552766j, -0.54307601-0.62108012j,\r\n       -0.23769241+2.91680425j])\r\n>>> cupy.linalg.solve(a, b)\r\narray([ 0.07332535+0.27916782j, -1.61809337-0.62108012j,\r\n       -0.23769241+2.91680425j])\r\n>>> numpy.linalg.solve(cupy.asnumpy(a), cupy.asnumpy(b))\r\narray([0.04332182-0.34957548j, 1.6260605 -0.2383031j ,\r\n       0.88753974+1.15498982j])\r\n>>> cupy.__version__\r\n'5.0.0b4'\r\n```\r\n\r\n#1524 would fix the issue.\n", "before_files": [{"content": "import numpy\nfrom numpy import linalg\nimport six\n\nimport cupy\nfrom cupy.core import core\nfrom cupy import cuda\nfrom cupy.cuda import cublas\nfrom cupy.cuda import device\nfrom cupy.linalg import decomposition\nfrom cupy.linalg import util\n\nif cuda.cusolver_enabled:\n    from cupy.cuda import cusolver\n\n\ndef solve(a, b):\n    \"\"\"Solves a linear matrix equation.\n\n    It computes the exact solution of ``x`` in ``ax = b``,\n    where ``a`` is a square and full rank matrix.\n\n    Args:\n        a (cupy.ndarray): The matrix with dimension ``(..., M, M)``.\n        b (cupy.ndarray): The matrix with dimension ``(...,M)`` or\n            ``(..., M, K)``.\n\n    Returns:\n        cupy.ndarray:\n            The matrix with dimension ``(..., M)`` or ``(..., M, K)``.\n\n    .. seealso:: :func:`numpy.linalg.solve`\n    \"\"\"\n    # NOTE: Since cusolver in CUDA 8.0 does not support gesv,\n    #       we manually solve a linear system with QR decomposition.\n    #       For details, please see the following:\n    #       https://docs.nvidia.com/cuda/cusolver/index.html#qr_examples\n    if not cuda.cusolver_enabled:\n        raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')\n\n    util._assert_cupy_array(a, b)\n    util._assert_nd_squareness(a)\n\n    if not ((a.ndim == b.ndim or a.ndim == b.ndim + 1) and\n            a.shape[:-1] == b.shape[:a.ndim - 1]):\n        raise ValueError(\n            'a must have (..., M, M) shape and b must have (..., M) '\n            'or (..., M, K)')\n\n    # Cast to float32 or float64\n    if a.dtype.char == 'f' or a.dtype.char == 'd':\n        dtype = a.dtype\n    else:\n        dtype = numpy.find_common_type((a.dtype.char, 'f'), ())\n\n    a = a.astype(dtype)\n    b = b.astype(dtype)\n    if a.ndim == 2:\n        return _solve(a, b)\n    x = cupy.empty_like(b)\n    shape = a.shape[:-2]\n    for i in six.moves.range(numpy.prod(shape)):\n        index = numpy.unravel_index(i, shape)\n        x[index] = _solve(a[index], b[index])\n    return x\n\n\ndef _solve(a, b):\n    a = cupy.asfortranarray(a)\n    b = cupy.asfortranarray(b)\n    dtype = a.dtype\n    m, k = (b.size, 1) if b.ndim == 1 else b.shape\n    cusolver_handle = device.get_cusolver_handle()\n    cublas_handle = device.get_cublas_handle()\n    dev_info = cupy.empty(1, dtype=numpy.int32)\n\n    if dtype == 'f':\n        geqrf = cusolver.sgeqrf\n        geqrf_bufferSize = cusolver.sgeqrf_bufferSize\n        ormqr = cusolver.sormqr\n        trsm = cublas.strsm\n    else:  # dtype == 'd'\n        geqrf = cusolver.dgeqrf\n        geqrf_bufferSize = cusolver.dgeqrf_bufferSize\n        ormqr = cusolver.dormqr\n        trsm = cublas.dtrsm\n\n    # 1. QR decomposition (A = Q * R)\n    buffersize = geqrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)\n    workspace = cupy.empty(buffersize, dtype=dtype)\n    tau = cupy.empty(m, dtype=dtype)\n    geqrf(\n        cusolver_handle, m, m, a.data.ptr, m,\n        tau.data.ptr, workspace.data.ptr, buffersize, dev_info.data.ptr)\n    _check_status(dev_info)\n    # 2. ormqr (Q^T * B)\n    ormqr(\n        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_OP_T,\n        m, k, m, a.data.ptr, m, tau.data.ptr, b.data.ptr, m,\n        workspace.data.ptr, buffersize, dev_info.data.ptr)\n    _check_status(dev_info)\n    # 3. trsm (X = R^{-1} * (Q^T * B))\n    trsm(\n        cublas_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_FILL_MODE_UPPER,\n        cublas.CUBLAS_OP_N, cublas.CUBLAS_DIAG_NON_UNIT,\n        m, k, 1, a.data.ptr, m, b.data.ptr, m)\n    return b\n\n\ndef _check_status(dev_info):\n    status = int(dev_info)\n    if status < 0:\n        raise linalg.LinAlgError(\n            'Parameter error (maybe caused by a bug in cupy.linalg?)')\n\n\ndef tensorsolve(a, b, axes=None):\n    \"\"\"Solves tensor equations denoted by ``ax = b``.\n\n    Suppose that ``b`` is equivalent to ``cupy.tensordot(a, x)``.\n    This function computes tensor ``x`` from ``a`` and ``b``.\n\n    Args:\n        a (cupy.ndarray): The tensor with ``len(shape) >= 1``\n        b (cupy.ndarray): The tensor with ``len(shape) >= 1``\n        axes (tuple of ints): Axes in ``a`` to reorder to the right\n            before inversion.\n\n    Returns:\n        cupy.ndarray:\n            The tensor with shape ``Q`` such that ``b.shape + Q == a.shape``.\n\n    .. seealso:: :func:`numpy.linalg.tensorsolve`\n    \"\"\"\n    if axes is not None:\n        allaxes = list(six.moves.range(a.ndim))\n        for k in axes:\n            allaxes.remove(k)\n            allaxes.insert(a.ndim, k)\n        a = a.transpose(allaxes)\n\n    oldshape = a.shape[-(a.ndim - b.ndim):]\n    prod = cupy.internal.prod(oldshape)\n\n    a = a.reshape(-1, prod)\n    b = b.ravel()\n    result = solve(a, b)\n    return result.reshape(oldshape)\n\n\n# TODO(okuta): Implement lstsq\n\n\ndef inv(a):\n    \"\"\"Computes the inverse of a matrix.\n\n    This function computes matrix ``a_inv`` from n-dimensional regular matrix\n    ``a`` such that ``dot(a, a_inv) == eye(n)``.\n\n    Args:\n        a (cupy.ndarray): The regular matrix\n\n    Returns:\n        cupy.ndarray: The inverse of a matrix.\n\n    .. seealso:: :func:`numpy.linalg.inv`\n    \"\"\"\n    if not cuda.cusolver_enabled:\n        raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')\n\n    # to prevent `a` to be overwritten\n    a = a.copy()\n\n    util._assert_cupy_array(a)\n    util._assert_rank2(a)\n    util._assert_nd_squareness(a)\n\n    if a.dtype.char == 'f' or a.dtype.char == 'd':\n        dtype = a.dtype.char\n    else:\n        dtype = numpy.find_common_type((a.dtype.char, 'f'), ()).char\n\n    cusolver_handle = device.get_cusolver_handle()\n    dev_info = cupy.empty(1, dtype=dtype)\n\n    ipiv = cupy.empty((a.shape[0], 1), dtype=dtype)\n\n    if dtype == 'f':\n        getrf = cusolver.sgetrf\n        getrf_bufferSize = cusolver.sgetrf_bufferSize\n        getrs = cusolver.sgetrs\n    else:  # dtype == 'd'\n        getrf = cusolver.dgetrf\n        getrf_bufferSize = cusolver.dgetrf_bufferSize\n        getrs = cusolver.dgetrs\n\n    m = a.shape[0]\n\n    buffersize = getrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)\n    workspace = cupy.empty(buffersize, dtype=dtype)\n\n    # LU factorization\n    getrf(cusolver_handle, m, m, a.data.ptr, m, workspace.data.ptr,\n          ipiv.data.ptr, dev_info.data.ptr)\n\n    b = cupy.eye(m, dtype=dtype)\n\n    # solve for the inverse\n    getrs(cusolver_handle, 0, m, m, a.data.ptr, m, ipiv.data.ptr, b.data.ptr,\n          m, dev_info.data.ptr)\n\n    return b\n\n\ndef pinv(a, rcond=1e-15):\n    \"\"\"Compute the Moore-Penrose pseudoinverse of a matrix.\n\n    It computes a pseudoinverse of a matrix ``a``, which is a generalization\n    of the inverse matrix with Singular Value Decomposition (SVD).\n    Note that it automatically removes small singular values for stability.\n\n    Args:\n        a (cupy.ndarray): The matrix with dimension ``(M, N)``\n        rcond (float): Cutoff parameter for small singular values.\n            For stability it computes the largest singular value denoted by\n            ``s``, and sets all singular values smaller than ``s`` to zero.\n\n    Returns:\n        cupy.ndarray: The pseudoinverse of ``a`` with dimension ``(N, M)``.\n\n    .. seealso:: :func:`numpy.linalg.pinv`\n    \"\"\"\n    u, s, vt = decomposition.svd(a, full_matrices=False)\n    cutoff = rcond * s.max()\n    s1 = 1 / s\n    s1[s <= cutoff] = 0\n    return core.dot(vt.T, s1[:, None] * u.T)\n\n\ndef tensorinv(a, ind=2):\n    \"\"\"Computes the inverse of a tensor.\n\n    This function computes tensor ``a_inv`` from tensor ``a`` such that\n    ``tensordot(a_inv, a, ind) == I``, where ``I`` denotes the identity tensor.\n\n    Args:\n        a (cupy.ndarray):\n            The tensor such that\n            ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.\n        ind (int):\n            The positive number used in ``axes`` option of ``tensordot``.\n\n    Returns:\n        cupy.ndarray:\n            The inverse of a tensor whose shape is equivalent to\n            ``a.shape[ind:] + a.shape[:ind]``.\n\n    .. seealso:: :func:`numpy.linalg.tensorinv`\n    \"\"\"\n    util._assert_cupy_array(a)\n\n    if ind <= 0:\n        raise ValueError('Invalid ind argument')\n    oldshape = a.shape\n    invshape = oldshape[ind:] + oldshape[:ind]\n    prod = cupy.internal.prod(oldshape[ind:])\n    a = a.reshape(prod, -1)\n    a_inv = inv(a)\n    return a_inv.reshape(*invshape)\n", "path": "cupy/linalg/solve.py"}], "after_files": [{"content": "import numpy\nfrom numpy import linalg\nimport six\n\nimport cupy\nfrom cupy.core import core\nfrom cupy import cuda\nfrom cupy.cuda import cublas\nfrom cupy.cuda import device\nfrom cupy.linalg import decomposition\nfrom cupy.linalg import util\n\nif cuda.cusolver_enabled:\n    from cupy.cuda import cusolver\n\n\ndef solve(a, b):\n    \"\"\"Solves a linear matrix equation.\n\n    It computes the exact solution of ``x`` in ``ax = b``,\n    where ``a`` is a square and full rank matrix.\n\n    Args:\n        a (cupy.ndarray): The matrix with dimension ``(..., M, M)``.\n        b (cupy.ndarray): The matrix with dimension ``(...,M)`` or\n            ``(..., M, K)``.\n\n    Returns:\n        cupy.ndarray:\n            The matrix with dimension ``(..., M)`` or ``(..., M, K)``.\n\n    .. seealso:: :func:`numpy.linalg.solve`\n    \"\"\"\n    # NOTE: Since cusolver in CUDA 8.0 does not support gesv,\n    #       we manually solve a linear system with QR decomposition.\n    #       For details, please see the following:\n    #       https://docs.nvidia.com/cuda/cusolver/index.html#qr_examples\n    if not cuda.cusolver_enabled:\n        raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')\n\n    util._assert_cupy_array(a, b)\n    util._assert_nd_squareness(a)\n\n    if not ((a.ndim == b.ndim or a.ndim == b.ndim + 1) and\n            a.shape[:-1] == b.shape[:a.ndim - 1]):\n        raise ValueError(\n            'a must have (..., M, M) shape and b must have (..., M) '\n            'or (..., M, K)')\n\n    # Cast to float32 or float64\n    if a.dtype.char == 'f' or a.dtype.char == 'd':\n        dtype = a.dtype\n    else:\n        dtype = numpy.find_common_type((a.dtype.char, 'f'), ())\n\n    cublas_handle = device.get_cublas_handle()\n    cusolver_handle = device.get_cusolver_handle()\n\n    a = a.astype(dtype)\n    b = b.astype(dtype)\n    if a.ndim == 2:\n        return _solve(a, b, cublas_handle, cusolver_handle)\n\n    x = cupy.empty_like(b)\n    shape = a.shape[:-2]\n    for i in six.moves.range(numpy.prod(shape)):\n        index = numpy.unravel_index(i, shape)\n        x[index] = _solve(a[index], b[index], cublas_handle, cusolver_handle)\n    return x\n\n\ndef _solve(a, b, cublas_handle, cusolver_handle):\n    a = cupy.asfortranarray(a)\n    b = cupy.asfortranarray(b)\n    dtype = a.dtype\n    m, k = (b.size, 1) if b.ndim == 1 else b.shape\n    dev_info = cupy.empty(1, dtype=numpy.int32)\n\n    if dtype == 'f':\n        geqrf = cusolver.sgeqrf\n        geqrf_bufferSize = cusolver.sgeqrf_bufferSize\n        ormqr = cusolver.sormqr\n        trans = cublas.CUBLAS_OP_T\n        trsm = cublas.strsm\n    elif dtype == 'd':\n        geqrf = cusolver.dgeqrf\n        geqrf_bufferSize = cusolver.dgeqrf_bufferSize\n        ormqr = cusolver.dormqr\n        trans = cublas.CUBLAS_OP_T\n        trsm = cublas.dtrsm\n    elif dtype == 'F':\n        geqrf = cusolver.cgeqrf\n        geqrf_bufferSize = cusolver.cgeqrf_bufferSize\n        ormqr = cusolver.cormqr\n        trans = cublas.CUBLAS_OP_C\n        trsm = cublas.ctrsm\n    elif dtype == 'D':\n        geqrf = cusolver.zgeqrf\n        geqrf_bufferSize = cusolver.zgeqrf_bufferSize\n        ormqr = cusolver.zormqr\n        trans = cublas.CUBLAS_OP_C\n        trsm = cublas.ztrsm\n    else:\n        raise NotImplementedError(dtype)\n\n    # 1. QR decomposition (A = Q * R)\n    buffersize = geqrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)\n    workspace = cupy.empty(buffersize, dtype=dtype)\n    tau = cupy.empty(m, dtype=dtype)\n    geqrf(\n        cusolver_handle, m, m, a.data.ptr, m,\n        tau.data.ptr, workspace.data.ptr, buffersize, dev_info.data.ptr)\n    _check_status(dev_info)\n    # 2. ormqr (Q^T * B)\n    ormqr(\n        cusolver_handle, cublas.CUBLAS_SIDE_LEFT, trans,\n        m, k, m, a.data.ptr, m, tau.data.ptr, b.data.ptr, m,\n        workspace.data.ptr, buffersize, dev_info.data.ptr)\n    _check_status(dev_info)\n    # 3. trsm (X = R^{-1} * (Q^T * B))\n    trsm(\n        cublas_handle, cublas.CUBLAS_SIDE_LEFT, cublas.CUBLAS_FILL_MODE_UPPER,\n        cublas.CUBLAS_OP_N, cublas.CUBLAS_DIAG_NON_UNIT,\n        m, k, 1, a.data.ptr, m, b.data.ptr, m)\n    return b\n\n\ndef _check_status(dev_info):\n    status = int(dev_info)\n    if status < 0:\n        raise linalg.LinAlgError(\n            'Parameter error (maybe caused by a bug in cupy.linalg?)')\n\n\ndef tensorsolve(a, b, axes=None):\n    \"\"\"Solves tensor equations denoted by ``ax = b``.\n\n    Suppose that ``b`` is equivalent to ``cupy.tensordot(a, x)``.\n    This function computes tensor ``x`` from ``a`` and ``b``.\n\n    Args:\n        a (cupy.ndarray): The tensor with ``len(shape) >= 1``\n        b (cupy.ndarray): The tensor with ``len(shape) >= 1``\n        axes (tuple of ints): Axes in ``a`` to reorder to the right\n            before inversion.\n\n    Returns:\n        cupy.ndarray:\n            The tensor with shape ``Q`` such that ``b.shape + Q == a.shape``.\n\n    .. seealso:: :func:`numpy.linalg.tensorsolve`\n    \"\"\"\n    if axes is not None:\n        allaxes = list(six.moves.range(a.ndim))\n        for k in axes:\n            allaxes.remove(k)\n            allaxes.insert(a.ndim, k)\n        a = a.transpose(allaxes)\n\n    oldshape = a.shape[-(a.ndim - b.ndim):]\n    prod = cupy.internal.prod(oldshape)\n\n    a = a.reshape(-1, prod)\n    b = b.ravel()\n    result = solve(a, b)\n    return result.reshape(oldshape)\n\n\n# TODO(okuta): Implement lstsq\n\n\ndef inv(a):\n    \"\"\"Computes the inverse of a matrix.\n\n    This function computes matrix ``a_inv`` from n-dimensional regular matrix\n    ``a`` such that ``dot(a, a_inv) == eye(n)``.\n\n    Args:\n        a (cupy.ndarray): The regular matrix\n\n    Returns:\n        cupy.ndarray: The inverse of a matrix.\n\n    .. seealso:: :func:`numpy.linalg.inv`\n    \"\"\"\n    if not cuda.cusolver_enabled:\n        raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')\n\n    # to prevent `a` to be overwritten\n    a = a.copy()\n\n    util._assert_cupy_array(a)\n    util._assert_rank2(a)\n    util._assert_nd_squareness(a)\n\n    if a.dtype.char == 'f' or a.dtype.char == 'd':\n        dtype = a.dtype.char\n    else:\n        dtype = numpy.find_common_type((a.dtype.char, 'f'), ()).char\n\n    cusolver_handle = device.get_cusolver_handle()\n    dev_info = cupy.empty(1, dtype=dtype)\n\n    ipiv = cupy.empty((a.shape[0], 1), dtype=dtype)\n\n    if dtype == 'f':\n        getrf = cusolver.sgetrf\n        getrf_bufferSize = cusolver.sgetrf_bufferSize\n        getrs = cusolver.sgetrs\n    else:  # dtype == 'd'\n        getrf = cusolver.dgetrf\n        getrf_bufferSize = cusolver.dgetrf_bufferSize\n        getrs = cusolver.dgetrs\n\n    m = a.shape[0]\n\n    buffersize = getrf_bufferSize(cusolver_handle, m, m, a.data.ptr, m)\n    workspace = cupy.empty(buffersize, dtype=dtype)\n\n    # LU factorization\n    getrf(cusolver_handle, m, m, a.data.ptr, m, workspace.data.ptr,\n          ipiv.data.ptr, dev_info.data.ptr)\n\n    b = cupy.eye(m, dtype=dtype)\n\n    # solve for the inverse\n    getrs(cusolver_handle, 0, m, m, a.data.ptr, m, ipiv.data.ptr, b.data.ptr,\n          m, dev_info.data.ptr)\n\n    return b\n\n\ndef pinv(a, rcond=1e-15):\n    \"\"\"Compute the Moore-Penrose pseudoinverse of a matrix.\n\n    It computes a pseudoinverse of a matrix ``a``, which is a generalization\n    of the inverse matrix with Singular Value Decomposition (SVD).\n    Note that it automatically removes small singular values for stability.\n\n    Args:\n        a (cupy.ndarray): The matrix with dimension ``(M, N)``\n        rcond (float): Cutoff parameter for small singular values.\n            For stability it computes the largest singular value denoted by\n            ``s``, and sets all singular values smaller than ``s`` to zero.\n\n    Returns:\n        cupy.ndarray: The pseudoinverse of ``a`` with dimension ``(N, M)``.\n\n    .. seealso:: :func:`numpy.linalg.pinv`\n    \"\"\"\n    u, s, vt = decomposition.svd(a, full_matrices=False)\n    cutoff = rcond * s.max()\n    s1 = 1 / s\n    s1[s <= cutoff] = 0\n    return core.dot(vt.T, s1[:, None] * u.T)\n\n\ndef tensorinv(a, ind=2):\n    \"\"\"Computes the inverse of a tensor.\n\n    This function computes tensor ``a_inv`` from tensor ``a`` such that\n    ``tensordot(a_inv, a, ind) == I``, where ``I`` denotes the identity tensor.\n\n    Args:\n        a (cupy.ndarray):\n            The tensor such that\n            ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.\n        ind (int):\n            The positive number used in ``axes`` option of ``tensordot``.\n\n    Returns:\n        cupy.ndarray:\n            The inverse of a tensor whose shape is equivalent to\n            ``a.shape[ind:] + a.shape[:ind]``.\n\n    .. seealso:: :func:`numpy.linalg.tensorinv`\n    \"\"\"\n    util._assert_cupy_array(a)\n\n    if ind <= 0:\n        raise ValueError('Invalid ind argument')\n    oldshape = a.shape\n    invshape = oldshape[ind:] + oldshape[:ind]\n    prod = cupy.internal.prod(oldshape[ind:])\n    a = a.reshape(prod, -1)\n    a_inv = inv(a)\n    return a_inv.reshape(*invshape)\n", "path": "cupy/linalg/solve.py"}]}