| 71 | // Moore-Penrose Pseudoinverse |
| 72 | template<typename T> |
| 73 | Array<T> pinverseSvd(const Array<T> &in, const double tol) { |
| 74 | in.eval(); |
| 75 | dim_t M = in.dims()[0]; |
| 76 | dim_t N = in.dims()[1]; |
| 77 | dim_t P = in.dims()[2]; |
| 78 | dim_t Q = in.dims()[3]; |
| 79 | |
| 80 | // Compute SVD |
| 81 | using Tr = typename dtype_traits<T>::base_type; |
| 82 | // Ideally, these initializations should use createEmptyArray(), but for |
| 83 | // some reason, linux-opencl-k80 will produce wrong results for large arrays |
| 84 | Array<T> u = createValueArray<T>(dim4(M, M, P, Q), scalar<T>(0)); |
| 85 | Array<T> vT = createValueArray<T>(dim4(N, N, P, Q), scalar<T>(0)); |
| 86 | Array<Tr> sVec = |
| 87 | createValueArray<Tr>(dim4(min(M, N), 1, P, Q), scalar<Tr>(0)); |
| 88 | for (dim_t j = 0; j < Q; ++j) { |
| 89 | for (dim_t i = 0; i < P; ++i) { |
| 90 | Array<T> inSlice = |
| 91 | getSubArray(in, false, 0, M - 1, 0, N - 1, i, i, j, j); |
| 92 | Array<Tr> sVecSlice = getSubArray( |
| 93 | sVec, false, 0, sVec.dims()[0] - 1, 0, 0, i, i, j, j); |
| 94 | Array<T> uSlice = getSubArray(u, false, 0, u.dims()[0] - 1, 0, |
| 95 | u.dims()[1] - 1, i, i, j, j); |
| 96 | Array<T> vTSlice = getSubArray(vT, false, 0, vT.dims()[0] - 1, 0, |
| 97 | vT.dims()[1] - 1, i, i, j, j); |
| 98 | svd<T, Tr>(sVecSlice, uSlice, vTSlice, inSlice); |
| 99 | } |
| 100 | } |
| 101 | |
| 102 | // Cast s back to original data type for matmul later |
| 103 | // (since svd() makes s' type the base type of T) |
| 104 | Array<T> sVecCast = cast<T, Tr>(sVec); |
| 105 | |
| 106 | Array<T> v = transpose(vT, true); |
| 107 | |
| 108 | // Build relative tolerance array |
| 109 | Array<Tr> sVecMax = reduce<af_max_t, Tr, Tr>(sVec, 0); |
| 110 | Array<T> sVecMaxCast = cast<T, Tr>(sVecMax); |
| 111 | double tolMulShape = tol * static_cast<double>(max(M, N)); |
| 112 | Array<T> tolMulShapeArr = |
| 113 | createValueArray<T>(sVecMaxCast.dims(), scalar<T>(tolMulShape)); |
| 114 | Array<T> relTol = |
| 115 | arithOp<T, af_mul_t>(tolMulShapeArr, sVecMaxCast, sVecMaxCast.dims()); |
| 116 | Array<T> relTolArr = tile<T>(relTol, dim4(sVecCast.dims()[0])); |
| 117 | |
| 118 | // Get reciprocal of sVec's non-zero values for s pinverse, except for |
| 119 | // very small non-zero values though (< relTol), in order to avoid very |
| 120 | // large reciprocals |
| 121 | Array<T> ones = createValueArray<T>(sVecCast.dims(), scalar<T>(1.)); |
| 122 | Array<T> sVecRecip = arithOp<T, af_div_t>(ones, sVecCast, sVecCast.dims()); |
| 123 | Array<char> cond = |
| 124 | logicOp<T, af_ge_t>(sVecCast, relTolArr, sVecCast.dims()); |
| 125 | Array<T> zeros = createValueArray<T>(sVecCast.dims(), scalar<T>(0.)); |
| 126 | sVecRecip = createSelectNode<T>(cond, sVecRecip, zeros, sVecRecip.dims()); |
| 127 | |
| 128 | // Make s vector into s pinverse array |
| 129 | Array<T> sVecRecipMod = modDims<T>( |
| 130 | sVecRecip, |