From 49cd34116f0ea17a12495a4d5284785566e8089a Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 25 Aug 2017 13:31:28 -0500 Subject: [PATCH 1/9] First pass at adding least-squares functionality for Python and Fortran (LAPACK) --- dagrt/builtins_python.py | 22 ++++++++++++ dagrt/codegen/fortran.py | 70 +++++++++++++++++++++++++++++++++++--- dagrt/function_registry.py | 36 ++++++++++++++++++++ 3 files changed, 124 insertions(+), 4 deletions(-) diff --git a/dagrt/builtins_python.py b/dagrt/builtins_python.py index 9f506cb..4dc64f3 100644 --- a/dagrt/builtins_python.py +++ b/dagrt/builtins_python.py @@ -103,6 +103,27 @@ def builtin_linear_solve(a, b, a_cols, b_cols): return res_mat.reshape(-1, order="F") +def builtin_lls(a, b, a_cols, a_rows, b_cols): + import numpy as np + if a_cols != np.floor(a_cols): + raise ValueError("lls() argument a_cols is not an integer") + if a_rows != np.floor(a_rows): + raise ValueError("lls() argument a_rows is not an integer") + if b_cols != np.floor(b_cols): + raise ValueError("lls() argument b_cols is not an integer") + a_cols = int(a_cols) + a_rows = int(a_rows) + b_cols = int(b_cols) + + a_mat = a.reshape(a_rows, a_cols, order="F") + b_mat = b.reshape(-1, b_cols, order="F") + + import numpy.linalg as la + res_mat = la.lstsq(a_mat, b_mat)[0] + + return res_mat.reshape(-1, order="F") + + def builtin_svd(a, a_cols): import numpy as np if a_cols != np.floor(a_cols): @@ -131,6 +152,7 @@ builtins = { "array": builtin_array, "matmul": builtin_matmul, "linear_solve": builtin_linear_solve, + "lls": builtin_lls, "svd": builtin_svd, "print": builtin_print, } diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index e8f6e85..b9c0931 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -1054,10 +1054,6 @@ class CodeGenerator(StructuredCodeGenerator): from .analysis import collect_ode_component_names_from_dag component_ids = collect_ode_component_names_from_dag(dag) - if not component_ids <= set(self.user_type_map): - raise RuntimeError("User type missing from user type map: %r" - % (component_ids - set(self.user_type_map))) - from dagrt.codegen.data import Scalar, UserType for comp_id in component_ids: self.sym_kind_table.set( @@ -2530,6 +2526,72 @@ builtin_linear_solve = CallCode(UTIL_MACROS + """ """) + +builtin_lls = CallCode(UTIL_MACROS + """ + <% + b_rows = declare_new("integer", "b_rows") + res_size = declare_new("integer", "res_size") + + %> + + ${b_rows} = int(${a_cols}) + ${res_size} = int(${a_cols}) + + <% + if a_kind != b_kind: + raise TypeError("lls requires both arguments " + "to have same kind") + + ltr = get_lapack_letter(a_kind) + + lls_temp = declare_new( + kind_to_fortran(a_kind)+", dimension(:), allocatable" + , "lls_temp") + work = declare_new( + kind_to_fortran(a_kind)+", dimension(:), allocatable" + , "work") + s = declare_new( + kind_to_fortran(a_kind)+", dimension(:), allocatable" + , "s") + info = declare_new("integer", "info") + lwork = declare_new("integer", "lwork") + rank = declare_new("integer", "rank") + rcond = declare_new("real", "rcond") + %> + + allocate(${lls_temp}(0:size(${a})-1)) + + ${lls_temp} = ${a} + ${rcond} = -1 + ${lwork} = 3*min(int(${a_rows}), int(${a_cols})) + max(2*min(int(${a_rows}), int(${a_cols})), max(int(${a_rows}), int(${a_cols})), int(${b_rows})) + + if (allocated(${result})) then + deallocate(${result}) + endif + + allocate(${result}(0:${res_size}-1)) + allocate(${s}(0:${res_size}-1)) + allocate(${work}(0:${lwork}-1)) + + ${result} = ${b} + + call ${ltr}gelss(int(${a_rows}), int(${a_cols}), & + int(${b_cols}), ${lls_temp}, int(${a_rows}), ${result}, & + int(${b_rows}), ${s}, ${rcond}, ${rank}, & + ${work}, ${lwork}, ${info}) + + if (${info}.ne.0) then + write(dagrt_stderr,*) & + 'gelss on ${a} failed with info=', ${info} + stop + endif + + deallocate(${lls_temp}) + deallocate(${s}) + deallocate(${work}) + + """) + builtin_print = CallCode(UTIL_MACROS + """ write(*,*) ${arg} """) diff --git a/dagrt/function_registry.py b/dagrt/function_registry.py index 7026b2d..6246633 100644 --- a/dagrt/function_registry.py +++ b/dagrt/function_registry.py @@ -353,6 +353,39 @@ class _LinearSolve(Function): return (Array(is_real_valued),) +class _LeastSquares(Function): + """``lls(a, b, a_cols, a_rows, b_cols)`` returns a 1D array containing the + matrix resulting from a linear least squares analysis + """ + + result_names = ("result",) + identifier = "lls" + arg_names = ("a", "b", "a_cols", "a_rows", "b_cols") + default_dict = {} + + def get_result_kinds(self, arg_kinds, check): + a_kind, b_kind, a_cols_kind, a_rows_kind, b_cols_kind = self.resolve_args(arg_kinds) + + if a_kind is None or b_kind is None: + raise UnableToInferKind( + "lls needs to know both arguments to infer result kind") + + if check and not isinstance(a_kind, Array): + raise TypeError("argument 'a' of 'lls' is not an array") + if check and not isinstance(b_kind, Array): + raise TypeError("argument 'a' of 'lls' is not an array") + if check and not isinstance(a_cols_kind, Scalar): + raise TypeError("argument 'a_cols' of 'lls' is not a scalar") + if check and not isinstance(a_rows_kind, Scalar): + raise TypeError("argument 'a_rows' of 'lls' is not a scalar") + if check and not isinstance(b_cols_kind, Scalar): + raise TypeError("argument 'b_cols' of 'lls' is not a scalar") + + is_real_valued = a_kind.is_real_valued and b_kind.is_real_valued + + return (Array(is_real_valued),) + + class _SVD(Function): """``linear_solve(a, a_cols)`` returns a 2D array ``u``, a 1D array ``sigma``, and a 2D array ``vt``, representing the (reduced) SVD of ``a``. @@ -425,6 +458,7 @@ def _make_bfr(): (_Array(), "self._builtin_array({args})"), (_MatMul(), "self._builtin_matmul({args})"), (_LinearSolve(), "self._builtin_linear_solve({args})"), + (_LeastSquares(), "self._builtin_lls({args})"), (_Print(), "self._builtin_print({args})"), (_SVD(), "self._builtin_svd({args})"), ]: @@ -450,6 +484,8 @@ def _make_bfr(): f.builtin_matmul) bfr = bfr.register_codegen(_LinearSolve.identifier, "fortran", f.builtin_linear_solve) + bfr = bfr.register_codegen(_LeastSquares.identifier, "fortran", + f.builtin_lls) bfr = bfr.register_codegen(_Print.identifier, "fortran", f.builtin_print) -- GitLab From a58ee503a724e6d98c7bc6fb3a7a6181b0a34f1a Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Wed, 30 Aug 2017 14:29:28 -0500 Subject: [PATCH 2/9] Fixes some bugs in DGELSS implementation --- dagrt/codegen/fortran.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index b9c0931..c2a7c8a 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -2534,7 +2534,7 @@ builtin_lls = CallCode(UTIL_MACROS + """ %> - ${b_rows} = int(${a_cols}) + ${check_matrix(b, b_cols, b_rows, "linear_solve")} ${res_size} = int(${a_cols}) <% @@ -2563,7 +2563,7 @@ builtin_lls = CallCode(UTIL_MACROS + """ ${lls_temp} = ${a} ${rcond} = -1 - ${lwork} = 3*min(int(${a_rows}), int(${a_cols})) + max(2*min(int(${a_rows}), int(${a_cols})), max(int(${a_rows}), int(${a_cols})), int(${b_rows})) + ${lwork} = 3*min(int(${a_rows}), int(${a_cols})) + max(2*min(int(${a_rows}), int(${a_cols})), max(int(${a_rows}), int(${a_cols})), 1) if (allocated(${result})) then deallocate(${result}) @@ -2573,11 +2573,11 @@ builtin_lls = CallCode(UTIL_MACROS + """ allocate(${s}(0:${res_size}-1)) allocate(${work}(0:${lwork}-1)) - ${result} = ${b} + ${result}(0:int(${b_rows})-1) = ${b} call ${ltr}gelss(int(${a_rows}), int(${a_cols}), & int(${b_cols}), ${lls_temp}, int(${a_rows}), ${result}, & - int(${b_rows}), ${s}, ${rcond}, ${rank}, & + int(${a_cols}), ${s}, ${rcond}, ${rank}, & ${work}, ${lwork}, ${info}) if (${info}.ne.0) then -- GitLab From a0675da4ebce8f063bc066aa03579637c36227ca Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 1 Sep 2017 16:35:13 -0500 Subject: [PATCH 3/9] Renaming stuff to match numpy convention --- dagrt/builtins_python.py | 10 +++---- dagrt/codegen/fortran.py | 14 +++++----- dagrt/function_registry.py | 20 +++++++------- test/test_codegen_fortran.py | 53 ++++++++++++++++++++++++++++++++++++ 4 files changed, 75 insertions(+), 22 deletions(-) diff --git a/dagrt/builtins_python.py b/dagrt/builtins_python.py index 4dc64f3..6b916b8 100644 --- a/dagrt/builtins_python.py +++ b/dagrt/builtins_python.py @@ -103,14 +103,14 @@ def builtin_linear_solve(a, b, a_cols, b_cols): return res_mat.reshape(-1, order="F") -def builtin_lls(a, b, a_cols, a_rows, b_cols): +def builtin_lstsq(a, b, a_cols, a_rows, b_cols): import numpy as np if a_cols != np.floor(a_cols): - raise ValueError("lls() argument a_cols is not an integer") + raise ValueError("lstsq() argument a_cols is not an integer") if a_rows != np.floor(a_rows): - raise ValueError("lls() argument a_rows is not an integer") + raise ValueError("lstsq() argument a_rows is not an integer") if b_cols != np.floor(b_cols): - raise ValueError("lls() argument b_cols is not an integer") + raise ValueError("lstsq() argument b_cols is not an integer") a_cols = int(a_cols) a_rows = int(a_rows) b_cols = int(b_cols) @@ -152,7 +152,7 @@ builtins = { "array": builtin_array, "matmul": builtin_matmul, "linear_solve": builtin_linear_solve, - "lls": builtin_lls, + "lstsq": builtin_lstsq, "svd": builtin_svd, "print": builtin_print, } diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index c2a7c8a..60813b6 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -2527,7 +2527,7 @@ builtin_linear_solve = CallCode(UTIL_MACROS + """ """) -builtin_lls = CallCode(UTIL_MACROS + """ +builtin_lstsq = CallCode(UTIL_MACROS + """ <% b_rows = declare_new("integer", "b_rows") res_size = declare_new("integer", "res_size") @@ -2539,12 +2539,12 @@ builtin_lls = CallCode(UTIL_MACROS + """ <% if a_kind != b_kind: - raise TypeError("lls requires both arguments " + raise TypeError("lstsq requires both arguments " "to have same kind") ltr = get_lapack_letter(a_kind) - lls_temp = declare_new( + lss_temp = declare_new( kind_to_fortran(a_kind)+", dimension(:), allocatable" , "lls_temp") work = declare_new( @@ -2559,9 +2559,9 @@ builtin_lls = CallCode(UTIL_MACROS + """ rcond = declare_new("real", "rcond") %> - allocate(${lls_temp}(0:size(${a})-1)) + allocate(${lss_temp}(0:size(${a})-1)) - ${lls_temp} = ${a} + ${lss_temp} = ${a} ${rcond} = -1 ${lwork} = 3*min(int(${a_rows}), int(${a_cols})) + max(2*min(int(${a_rows}), int(${a_cols})), max(int(${a_rows}), int(${a_cols})), 1) @@ -2576,7 +2576,7 @@ builtin_lls = CallCode(UTIL_MACROS + """ ${result}(0:int(${b_rows})-1) = ${b} call ${ltr}gelss(int(${a_rows}), int(${a_cols}), & - int(${b_cols}), ${lls_temp}, int(${a_rows}), ${result}, & + int(${b_cols}), ${lss_temp}, int(${a_rows}), ${result}, & int(${a_cols}), ${s}, ${rcond}, ${rank}, & ${work}, ${lwork}, ${info}) @@ -2586,7 +2586,7 @@ builtin_lls = CallCode(UTIL_MACROS + """ stop endif - deallocate(${lls_temp}) + deallocate(${lss_temp}) deallocate(${s}) deallocate(${work}) diff --git a/dagrt/function_registry.py b/dagrt/function_registry.py index 6246633..89360af 100644 --- a/dagrt/function_registry.py +++ b/dagrt/function_registry.py @@ -354,12 +354,12 @@ class _LinearSolve(Function): class _LeastSquares(Function): - """``lls(a, b, a_cols, a_rows, b_cols)`` returns a 1D array containing the + """``lstsq(a, b, a_cols, a_rows, b_cols)`` returns a 1D array containing the matrix resulting from a linear least squares analysis """ result_names = ("result",) - identifier = "lls" + identifier = "lstsq" arg_names = ("a", "b", "a_cols", "a_rows", "b_cols") default_dict = {} @@ -368,18 +368,18 @@ class _LeastSquares(Function): if a_kind is None or b_kind is None: raise UnableToInferKind( - "lls needs to know both arguments to infer result kind") + "lstsq needs to know both arguments to infer result kind") if check and not isinstance(a_kind, Array): - raise TypeError("argument 'a' of 'lls' is not an array") + raise TypeError("argument 'a' of 'lstsq' is not an array") if check and not isinstance(b_kind, Array): - raise TypeError("argument 'a' of 'lls' is not an array") + raise TypeError("argument 'a' of 'lstsq' is not an array") if check and not isinstance(a_cols_kind, Scalar): - raise TypeError("argument 'a_cols' of 'lls' is not a scalar") + raise TypeError("argument 'a_cols' of 'lstsq' is not a scalar") if check and not isinstance(a_rows_kind, Scalar): - raise TypeError("argument 'a_rows' of 'lls' is not a scalar") + raise TypeError("argument 'a_rows' of 'lstsq' is not a scalar") if check and not isinstance(b_cols_kind, Scalar): - raise TypeError("argument 'b_cols' of 'lls' is not a scalar") + raise TypeError("argument 'b_cols' of 'lstsq' is not a scalar") is_real_valued = a_kind.is_real_valued and b_kind.is_real_valued @@ -458,7 +458,7 @@ def _make_bfr(): (_Array(), "self._builtin_array({args})"), (_MatMul(), "self._builtin_matmul({args})"), (_LinearSolve(), "self._builtin_linear_solve({args})"), - (_LeastSquares(), "self._builtin_lls({args})"), + (_LeastSquares(), "self._builtin_lstsq({args})"), (_Print(), "self._builtin_print({args})"), (_SVD(), "self._builtin_svd({args})"), ]: @@ -485,7 +485,7 @@ def _make_bfr(): bfr = bfr.register_codegen(_LinearSolve.identifier, "fortran", f.builtin_linear_solve) bfr = bfr.register_codegen(_LeastSquares.identifier, "fortran", - f.builtin_lls) + f.builtin_lstsq) bfr = bfr.register_codegen(_Print.identifier, "fortran", f.builtin_print) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index 919f188..505a032 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -121,6 +121,59 @@ def test_arrays_and_linalg(): fortran_options=["-llapack", "-lblas"]) +def test_least_squares(): + from dagrt.function_registry import base_function_registry as freg + + with CodeBuilder(label="primary") as cb: + cb("n", "4") + cb("nodes", "`array`(n)") + cb("vdm", "`array`(n*n)") + cb("identity", "`array`(n*n)") + cb.fence() + + cb("nodes[i]", "i/n", + loops=[("i", 0, "n")]) + cb("identity[i]", "0", + loops=[("i", 0, "n*n")]) + cb.fence() + + cb("identity[i*n + i]", "1", + loops=[("i", 0, "n")]) + cb("vdm[j*n + i]", "nodes[i]**j", + loops=[("i", 0, "n"), ("j", 0, "n")]) + + cb.fence() + + cb("vdm_inverse", "`lls`(vdm, identity, n, n)") + cb("myarray", "`matmul`(vdm, vdm_inverse, n, n)") + + cb("myzero", "myarray - identity") + cb((), "`print`(myzero)") + with cb.if_("`norm_2`(myzero) > 10**(-8)"): + cb.raise_(MatrixInversionFailure) + + code = DAGCode.create_with_steady_state( + cb.state_dependencies, cb.instructions) + + codegen = f.CodeGenerator( + 'arrays', + function_registry=freg, + user_type_map={}, + emit_instrumentation=True, + timing_function="second") + + code_str = codegen(code) + if 0: + with open("arrays.f90", "wt") as outf: + outf.write(code_str) + + run_fortran([ + ("arrays.f90", code_str), + ("test_arrays_and_linalg.f90", read_file("test_arrays_and_linalg.f90")), + ], + fortran_options=["-llapack", "-lblas"]) + + if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) -- GitLab From 37f12107640f0598e6afa9ba39132d595b0e4fbc Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 19 Sep 2017 10:32:27 -0500 Subject: [PATCH 4/9] Work towards SVD builtin for Fortran --- dagrt/builtins_python.py | 14 ++++++ dagrt/codegen/fortran.py | 97 +++++++++++++++++++++++++++++++++++++- dagrt/function_registry.py | 38 +++++++++++++-- 3 files changed, 145 insertions(+), 4 deletions(-) diff --git a/dagrt/builtins_python.py b/dagrt/builtins_python.py index 6b916b8..5570eb1 100644 --- a/dagrt/builtins_python.py +++ b/dagrt/builtins_python.py @@ -85,6 +85,19 @@ def builtin_matmul(a, b, a_cols, b_cols): return res_mat.reshape(-1, order="F") +def builtin_transpose(a, a_cols): + import numpy as np + if a_cols != np.floor(a_cols): + raise ValueError("transpose() argument a_cols is not an integer") + a_cols = int(a_cols) + + a_mat = a.reshape(-1, a_cols, order="F") + + res_mat = np.transpose(a_mat) + + return res_mat.reshape(-1, order="F") + + def builtin_linear_solve(a, b, a_cols, b_cols): import numpy as np if a_cols != np.floor(a_cols): @@ -151,6 +164,7 @@ builtins = { "dot_product": builtin_dot_product, "array": builtin_array, "matmul": builtin_matmul, + "transpose": builtin_transpose, "linear_solve": builtin_linear_solve, "lstsq": builtin_lstsq, "svd": builtin_svd, diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index 60813b6..1cfb42d 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -2400,7 +2400,7 @@ UTIL_MACROS = """ if (${rows_var} * int(${cols_var}) .ne. size(${mat_array})) then write(dagrt_stderr,*) & 'size of argument ' // & - '${mat_array}' // & + '${mat_array} ' // & 'to ${func_name} ' // & 'not divisible by ' // & '${cols_var}' @@ -2463,6 +2463,31 @@ builtin_matmul = CallCode(UTIL_MACROS + """ (/${res_size}/)) """) + +builtin_transpose = CallCode(UTIL_MACROS + """ + <% + a_rows = declare_new("integer", "a_rows") + res_size = declare_new("integer", "res_size") + %> + + ${check_matrix(a, a_cols, a_rows, "transpose")} + + ${a_rows} = size(${a}) / int(${a_cols}) + ${res_size} = ${a_rows} * int(${a_cols}) + + if (allocated(${result})) then + deallocate(${result}) + endif + + allocate(${result}(0:${res_size}-1)) + + ${result} = reshape( & + transpose( & + reshape(${a}, (/${a_rows}, int(${a_cols})/))), & + (/${res_size}/)) + """) + + builtin_linear_solve = CallCode(UTIL_MACROS + """ <% a_rows = declare_new("integer", "a_rows") @@ -2527,6 +2552,69 @@ builtin_linear_solve = CallCode(UTIL_MACROS + """ """) +builtin_svd = CallCode(UTIL_MACROS + """ + <% + sigma_size = declare_new("integer", "res_size") + a_rows = declare_new("integer", "a_rows") + + %> + + ${check_matrix(a, a_cols, a_rows, "svd")} + ${sigma_size} = min(int(${a_cols}),int(${a_rows})) + + <% + ltr = get_lapack_letter(a_kind) + + a_temp = declare_new( + kind_to_fortran(a_kind)+", dimension(:), allocatable" + , "a_temp") + work = declare_new( + kind_to_fortran(a_kind)+", dimension(:), allocatable" + , "work") + info = declare_new("integer", "info") + lwork = declare_new("integer", "lwork") + lda = declare_new("integer", "lda") + ldu = declare_new("integer", "ldu") + ldvt = declare_new("integer", "ldvt") + jobu = declare_new("character*1", "jobu") + jobvt = declare_new("character*1", "jobvt") + %> + + allocate(${a_temp}(0:size(${a})-1)) + ${jobu} = "S" + ${jobvt} = "S" + ${lda} = max(1,int(${a_rows})) + ${ldu} = int(${a_rows}) + ${ldvt} = min(int(${a_rows}),int(${a_rows})) + + ${a_temp} = ${a} + ${lwork} = max(1,3*min(int(${a_rows}), int(${a_cols})) + max(int(${a_rows}), int(${a_cols})), 5*min(int(${a_rows}), int(${a_cols}))) + + if (allocated(${sigma})) then + deallocate(${sigma}) + endif + + allocate(${sigma}(0:${sigma_size}-1)) + allocate(${work}(0:${lwork}-1)) + allocate(${u}(0:int(${a_rows}*${a_rows})-1)) + allocate(${vt}(0:int(${a_rows}*${a_cols})-1)) + + call ${ltr}gesvd(${jobu}, ${jobvt}, & + int(${a_rows}), int(${a_cols}), ${a_temp}, ${lda}, ${sigma}, & + ${u}, ${ldu}, ${vt}, ${ldvt}, ${work}, ${lwork}, ${info}) + + if (${info}.ne.0) then + write(dagrt_stderr,*) & + 'gesvd on ${a} failed with info=', ${info} + stop + endif + + deallocate(${a_temp}) + deallocate(${work}) + + """) + + builtin_lstsq = CallCode(UTIL_MACROS + """ <% b_rows = declare_new("integer", "b_rows") @@ -2557,9 +2645,11 @@ builtin_lstsq = CallCode(UTIL_MACROS + """ lwork = declare_new("integer", "lwork") rank = declare_new("integer", "rank") rcond = declare_new("real", "rcond") + trans = declare_new("character*1", "trans") %> allocate(${lss_temp}(0:size(${a})-1)) + ${trans} = "N" ${lss_temp} = ${a} ${rcond} = -1 @@ -2574,12 +2664,17 @@ builtin_lstsq = CallCode(UTIL_MACROS + """ allocate(${work}(0:${lwork}-1)) ${result}(0:int(${b_rows})-1) = ${b} + ${result}(int(${b_rows})) = 0 call ${ltr}gelss(int(${a_rows}), int(${a_cols}), & int(${b_cols}), ${lss_temp}, int(${a_rows}), ${result}, & int(${a_cols}), ${s}, ${rcond}, ${rank}, & ${work}, ${lwork}, ${info}) + !call ${ltr}gels(${trans}, int(${a_rows}), int(${a_cols}), & + ! int(${b_cols}), ${lss_temp}, int(${a_rows}), ${result}, & + ! int(${a_cols}), ${work}, ${lwork}, ${info}) + if (${info}.ne.0) then write(dagrt_stderr,*) & 'gelss on ${a} failed with info=', ${info} diff --git a/dagrt/function_registry.py b/dagrt/function_registry.py index 89360af..5228b3b 100644 --- a/dagrt/function_registry.py +++ b/dagrt/function_registry.py @@ -321,6 +321,33 @@ class _MatMul(Function): return (Array(is_real_valued),) +class _Transpose(Function): + """``transpose(a, a_cols)`` returns a 1D array containing the + matrix resulting from transposing the array a + """ + + result_names = ("result",) + identifier = "transpose" + arg_names = ("a", "a_cols") + default_dict = {} + + def get_result_kinds(self, arg_kinds, check): + a_kind, a_cols_kind = self.resolve_args(arg_kinds) + + if a_kind is None: + raise UnableToInferKind( + "transpose needs to know both arguments to infer result kind") + + if check and not isinstance(a_kind, Array): + raise TypeError("argument 'a' of 'transpose' is not an array") + if check and not isinstance(a_cols_kind, Scalar): + raise TypeError("argument 'a_cols' of 'transpose' is not a scalar") + + is_real_valued = a_kind.is_real_valued + + return (Array(is_real_valued),) + + class _LinearSolve(Function): """``linear_solve(a, b, a_cols, b_cols)`` returns a 1D array containing the matrix resulting from multiplying the matrix inverse of a by b (both interpreted @@ -387,7 +414,7 @@ class _LeastSquares(Function): class _SVD(Function): - """``linear_solve(a, a_cols)`` returns a 2D array ``u``, a 1D array ``sigma``, and + """``SVD(a, a_cols)`` returns a 2D array ``u``, a 1D array ``sigma``, and a 2D array ``vt``, representing the (reduced) SVD of ``a``. """ @@ -404,9 +431,9 @@ class _SVD(Function): "svd needs to know its argument to infer result kind") if check and not isinstance(a_kind, Array): - raise TypeError("argument 'a' of 'linear_solve' is not an array") + raise TypeError("argument 'a' of 'svd' is not an array") if check and not isinstance(a_cols_kind, Scalar): - raise TypeError("argument 'a_cols' of 'linear_solve' is not a scalar") + raise TypeError("argument 'a_cols' of 'svd' is not a scalar") is_real_valued = a_kind.is_real_valued @@ -457,6 +484,7 @@ def _make_bfr(): (_IsNaN(), "{numpy}.isnan({args})"), (_Array(), "self._builtin_array({args})"), (_MatMul(), "self._builtin_matmul({args})"), + (_Transpose(), "self._builtin_transpose({args})"), (_LinearSolve(), "self._builtin_linear_solve({args})"), (_LeastSquares(), "self._builtin_lstsq({args})"), (_Print(), "self._builtin_print({args})"), @@ -482,10 +510,14 @@ def _make_bfr(): f.builtin_array) bfr = bfr.register_codegen(_MatMul.identifier, "fortran", f.builtin_matmul) + bfr = bfr.register_codegen(_Transpose.identifier, "fortran", + f.builtin_transpose) bfr = bfr.register_codegen(_LinearSolve.identifier, "fortran", f.builtin_linear_solve) bfr = bfr.register_codegen(_LeastSquares.identifier, "fortran", f.builtin_lstsq) + bfr = bfr.register_codegen(_SVD.identifier, "fortran", + f.builtin_svd) bfr = bfr.register_codegen(_Print.identifier, "fortran", f.builtin_print) -- GitLab From cb177255d0e0061bcf938794540b5e45e14d3641 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 19 Sep 2017 11:05:52 -0500 Subject: [PATCH 5/9] Placate flake8 --- dagrt/codegen/fortran.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index 1cfb42d..1227f80 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -47,6 +47,7 @@ def pad_fortran(line, width): line += '&' return line + wrap_line = partial(wrap_line_base, pad_func=pad_fortran) @@ -2588,7 +2589,11 @@ builtin_svd = CallCode(UTIL_MACROS + """ ${ldvt} = min(int(${a_rows}),int(${a_rows})) ${a_temp} = ${a} - ${lwork} = max(1,3*min(int(${a_rows}), int(${a_cols})) + max(int(${a_rows}), int(${a_cols})), 5*min(int(${a_rows}), int(${a_cols}))) + ${lwork} = max(1, & + 3*min(int(${a_rows}), & + int(${a_cols})) + max(int(${a_rows}), & + int(${a_cols})), & + 5*min(int(${a_rows}), int(${a_cols}))) if (allocated(${sigma})) then deallocate(${sigma}) @@ -2653,7 +2658,12 @@ builtin_lstsq = CallCode(UTIL_MACROS + """ ${lss_temp} = ${a} ${rcond} = -1 - ${lwork} = 3*min(int(${a_rows}), int(${a_cols})) + max(2*min(int(${a_rows}), int(${a_cols})), max(int(${a_rows}), int(${a_cols})), 1) + ${lwork} = 3*min( & + int(${a_rows}), & + int(${a_cols})) + max(2*min(int(${a_rows}), & + int(${a_cols})), & + max(int(${a_rows}), & + int(${a_cols})), 1) if (allocated(${result})) then deallocate(${result}) -- GitLab From bd5ed453b4961d7e8bf64009d3dda23e6b65a0c1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 19 Sep 2017 11:10:36 -0500 Subject: [PATCH 6/9] Fix codegen for function calls with multiple return values --- dagrt/codegen/fortran.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index 1227f80..6922a8e 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -2151,7 +2151,7 @@ class CodeGenerator(StructuredCodeGenerator): inst.as_expression()) assignee_fortran_names = [ - self.name_manager[assignee_sym] for a in inst.assignees] + self.name_manager[assignee_sym] for assignee_sym in inst.assignees] function = self.function_registry[inst.function_id] -- GitLab From 1381cf0b61a34f1cf5541862f383eac565c5c2f2 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 10:54:00 -0500 Subject: [PATCH 7/9] Remove outdated and ineffectual least squares test --- test/test_codegen_fortran.py | 53 ------------------------------------ 1 file changed, 53 deletions(-) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index 505a032..919f188 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -121,59 +121,6 @@ def test_arrays_and_linalg(): fortran_options=["-llapack", "-lblas"]) -def test_least_squares(): - from dagrt.function_registry import base_function_registry as freg - - with CodeBuilder(label="primary") as cb: - cb("n", "4") - cb("nodes", "`array`(n)") - cb("vdm", "`array`(n*n)") - cb("identity", "`array`(n*n)") - cb.fence() - - cb("nodes[i]", "i/n", - loops=[("i", 0, "n")]) - cb("identity[i]", "0", - loops=[("i", 0, "n*n")]) - cb.fence() - - cb("identity[i*n + i]", "1", - loops=[("i", 0, "n")]) - cb("vdm[j*n + i]", "nodes[i]**j", - loops=[("i", 0, "n"), ("j", 0, "n")]) - - cb.fence() - - cb("vdm_inverse", "`lls`(vdm, identity, n, n)") - cb("myarray", "`matmul`(vdm, vdm_inverse, n, n)") - - cb("myzero", "myarray - identity") - cb((), "`print`(myzero)") - with cb.if_("`norm_2`(myzero) > 10**(-8)"): - cb.raise_(MatrixInversionFailure) - - code = DAGCode.create_with_steady_state( - cb.state_dependencies, cb.instructions) - - codegen = f.CodeGenerator( - 'arrays', - function_registry=freg, - user_type_map={}, - emit_instrumentation=True, - timing_function="second") - - code_str = codegen(code) - if 0: - with open("arrays.f90", "wt") as outf: - outf.write(code_str) - - run_fortran([ - ("arrays.f90", code_str), - ("test_arrays_and_linalg.f90", read_file("test_arrays_and_linalg.f90")), - ], - fortran_options=["-llapack", "-lblas"]) - - if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) -- GitLab From c3b9b363da775477f678ca1a4d511998701f2e36 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 12:46:11 -0500 Subject: [PATCH 8/9] Removes outdated lstsq stuff --- dagrt/builtins_python.py | 22 ----------- dagrt/codegen/fortran.py | 77 -------------------------------------- dagrt/function_registry.py | 37 +----------------- 3 files changed, 1 insertion(+), 135 deletions(-) diff --git a/dagrt/builtins_python.py b/dagrt/builtins_python.py index 5570eb1..623776c 100644 --- a/dagrt/builtins_python.py +++ b/dagrt/builtins_python.py @@ -116,27 +116,6 @@ def builtin_linear_solve(a, b, a_cols, b_cols): return res_mat.reshape(-1, order="F") -def builtin_lstsq(a, b, a_cols, a_rows, b_cols): - import numpy as np - if a_cols != np.floor(a_cols): - raise ValueError("lstsq() argument a_cols is not an integer") - if a_rows != np.floor(a_rows): - raise ValueError("lstsq() argument a_rows is not an integer") - if b_cols != np.floor(b_cols): - raise ValueError("lstsq() argument b_cols is not an integer") - a_cols = int(a_cols) - a_rows = int(a_rows) - b_cols = int(b_cols) - - a_mat = a.reshape(a_rows, a_cols, order="F") - b_mat = b.reshape(-1, b_cols, order="F") - - import numpy.linalg as la - res_mat = la.lstsq(a_mat, b_mat)[0] - - return res_mat.reshape(-1, order="F") - - def builtin_svd(a, a_cols): import numpy as np if a_cols != np.floor(a_cols): @@ -166,7 +145,6 @@ builtins = { "matmul": builtin_matmul, "transpose": builtin_transpose, "linear_solve": builtin_linear_solve, - "lstsq": builtin_lstsq, "svd": builtin_svd, "print": builtin_print, } diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index 6922a8e..c25df4d 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -2620,83 +2620,6 @@ builtin_svd = CallCode(UTIL_MACROS + """ """) -builtin_lstsq = CallCode(UTIL_MACROS + """ - <% - b_rows = declare_new("integer", "b_rows") - res_size = declare_new("integer", "res_size") - - %> - - ${check_matrix(b, b_cols, b_rows, "linear_solve")} - ${res_size} = int(${a_cols}) - - <% - if a_kind != b_kind: - raise TypeError("lstsq requires both arguments " - "to have same kind") - - ltr = get_lapack_letter(a_kind) - - lss_temp = declare_new( - kind_to_fortran(a_kind)+", dimension(:), allocatable" - , "lls_temp") - work = declare_new( - kind_to_fortran(a_kind)+", dimension(:), allocatable" - , "work") - s = declare_new( - kind_to_fortran(a_kind)+", dimension(:), allocatable" - , "s") - info = declare_new("integer", "info") - lwork = declare_new("integer", "lwork") - rank = declare_new("integer", "rank") - rcond = declare_new("real", "rcond") - trans = declare_new("character*1", "trans") - %> - - allocate(${lss_temp}(0:size(${a})-1)) - ${trans} = "N" - - ${lss_temp} = ${a} - ${rcond} = -1 - ${lwork} = 3*min( & - int(${a_rows}), & - int(${a_cols})) + max(2*min(int(${a_rows}), & - int(${a_cols})), & - max(int(${a_rows}), & - int(${a_cols})), 1) - - if (allocated(${result})) then - deallocate(${result}) - endif - - allocate(${result}(0:${res_size}-1)) - allocate(${s}(0:${res_size}-1)) - allocate(${work}(0:${lwork}-1)) - - ${result}(0:int(${b_rows})-1) = ${b} - ${result}(int(${b_rows})) = 0 - - call ${ltr}gelss(int(${a_rows}), int(${a_cols}), & - int(${b_cols}), ${lss_temp}, int(${a_rows}), ${result}, & - int(${a_cols}), ${s}, ${rcond}, ${rank}, & - ${work}, ${lwork}, ${info}) - - !call ${ltr}gels(${trans}, int(${a_rows}), int(${a_cols}), & - ! int(${b_cols}), ${lss_temp}, int(${a_rows}), ${result}, & - ! int(${a_cols}), ${work}, ${lwork}, ${info}) - - if (${info}.ne.0) then - write(dagrt_stderr,*) & - 'gelss on ${a} failed with info=', ${info} - stop - endif - - deallocate(${lss_temp}) - deallocate(${s}) - deallocate(${work}) - - """) - builtin_print = CallCode(UTIL_MACROS + """ write(*,*) ${arg} """) diff --git a/dagrt/function_registry.py b/dagrt/function_registry.py index 5228b3b..5082151 100644 --- a/dagrt/function_registry.py +++ b/dagrt/function_registry.py @@ -380,39 +380,6 @@ class _LinearSolve(Function): return (Array(is_real_valued),) -class _LeastSquares(Function): - """``lstsq(a, b, a_cols, a_rows, b_cols)`` returns a 1D array containing the - matrix resulting from a linear least squares analysis - """ - - result_names = ("result",) - identifier = "lstsq" - arg_names = ("a", "b", "a_cols", "a_rows", "b_cols") - default_dict = {} - - def get_result_kinds(self, arg_kinds, check): - a_kind, b_kind, a_cols_kind, a_rows_kind, b_cols_kind = self.resolve_args(arg_kinds) - - if a_kind is None or b_kind is None: - raise UnableToInferKind( - "lstsq needs to know both arguments to infer result kind") - - if check and not isinstance(a_kind, Array): - raise TypeError("argument 'a' of 'lstsq' is not an array") - if check and not isinstance(b_kind, Array): - raise TypeError("argument 'a' of 'lstsq' is not an array") - if check and not isinstance(a_cols_kind, Scalar): - raise TypeError("argument 'a_cols' of 'lstsq' is not a scalar") - if check and not isinstance(a_rows_kind, Scalar): - raise TypeError("argument 'a_rows' of 'lstsq' is not a scalar") - if check and not isinstance(b_cols_kind, Scalar): - raise TypeError("argument 'b_cols' of 'lstsq' is not a scalar") - - is_real_valued = a_kind.is_real_valued and b_kind.is_real_valued - - return (Array(is_real_valued),) - - class _SVD(Function): """``SVD(a, a_cols)`` returns a 2D array ``u``, a 1D array ``sigma``, and a 2D array ``vt``, representing the (reduced) SVD of ``a``. @@ -486,7 +453,6 @@ def _make_bfr(): (_MatMul(), "self._builtin_matmul({args})"), (_Transpose(), "self._builtin_transpose({args})"), (_LinearSolve(), "self._builtin_linear_solve({args})"), - (_LeastSquares(), "self._builtin_lstsq({args})"), (_Print(), "self._builtin_print({args})"), (_SVD(), "self._builtin_svd({args})"), ]: @@ -514,8 +480,6 @@ def _make_bfr(): f.builtin_transpose) bfr = bfr.register_codegen(_LinearSolve.identifier, "fortran", f.builtin_linear_solve) - bfr = bfr.register_codegen(_LeastSquares.identifier, "fortran", - f.builtin_lstsq) bfr = bfr.register_codegen(_SVD.identifier, "fortran", f.builtin_svd) bfr = bfr.register_codegen(_Print.identifier, "fortran", @@ -523,6 +487,7 @@ def _make_bfr(): return bfr + base_function_registry = _make_bfr() # }}} -- GitLab From 24b064863d847f0a656af46ab7b319d9258823d6 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 15:59:46 -0500 Subject: [PATCH 9/9] Fixes errantly deleted if statement --- dagrt/codegen/fortran.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dagrt/codegen/fortran.py b/dagrt/codegen/fortran.py index c25df4d..fe4f3b7 100644 --- a/dagrt/codegen/fortran.py +++ b/dagrt/codegen/fortran.py @@ -1055,6 +1055,10 @@ class CodeGenerator(StructuredCodeGenerator): from .analysis import collect_ode_component_names_from_dag component_ids = collect_ode_component_names_from_dag(dag) + if not component_ids <= set(self.user_type_map): + raise RuntimeError("User type missing from user type map: %r" + % (component_ids - set(self.user_type_map))) + from dagrt.codegen.data import Scalar, UserType for comp_id in component_ids: self.sym_kind_table.set( -- GitLab