Release 5.0.0-alpha.

Merge pull request #26202 from gursimarsingh:improved_tracker_samples
Improved Tracker Samples #26202 Relates to #25006 This sample has been rewritten to track a selected target in a video or camera stream. It combines VIT tracker, Nano tracker and Dasiamrpn tracker into one tracker sample ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [x] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake
2026-07-25 21:33:04 +04:00 · 2024-12-05 11:54:05 +03:00 · 2024-12-05 11:50:03 +03:00 · 2024-12-05 11:34:49 +03:00 · 2024-12-05 09:10:41 +03:00 · 2024-12-05 09:03:51 +03:00
3546 changed files with 202533 additions and 1198943 deletions
@@ -1,14 +0,0 @@
-name: 4.x
-
-on:
-  schedule:
-    - cron: '0 3 * * *'
-  workflow_dispatch:
-
-jobs:
-  CodeQL:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-CodeQL.yaml@main
-    with:
-      target_branch: '4.x'
-      workflow_branch: main
-
@@ -1,55 +0,0 @@
-name: PR:4.x
-
-on:
-  pull_request:
-    branches:
-      - 4.x
-
-jobs:
-
-  Linux:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-Linux.yaml@main
-    with:
-      workflow_branch: main
-
-  Windows:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-Windows.yaml@main
-    with:
-      workflow_branch: main
-
-  Ubuntu2404-ARM64:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-ARM64.yaml@main
-
-  Ubuntu2404-ARM64-Debug:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-ARM64-Debug.yaml@main
-
-  Ubuntu2004-x64-OpenVINO:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-U20-OpenVINO.yaml@main
-
-  Ubuntu2004-x64-CUDA:
-    if: "${{ contains(github.event.pull_request.labels.*.name, 'category: dnn') }} || ${{ contains(github.event.pull_request.labels.*.name, 'category: dnn (onnx)') }}"
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-U20-Cuda.yaml@main
-
-  macOS-ARM64:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-macOS-ARM64.yaml@main
-
-  macOS-x64:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-macOS-x86_64.yaml@main
-
-  macOS-ARM64-Vulkan:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-macOS-ARM64-Vulkan.yaml@main
-
-  iOS:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-iOS.yaml@main
-
-  Android-SDK:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-4.x-Android-SDK.yaml@main
-
-  TIM-VX:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-timvx-backend-tests-4.x.yml@main
-
-  docs:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-docs.yaml@main
-
-  Linux-RISC-V-Clang:
-    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-4.x-RISCV.yaml@main
@@ -0,0 +1,67 @@
+name: PR:5.x
+
+on:
+  pull_request:
+    branches:
+      - 5.x
+
+jobs:
+  Ubuntu2004-ARM64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-ARM64.yaml@main
+
+  Ubuntu2004-ARM64-Debug:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-ARM64-Debug.yaml@main
+
+  Ubuntu2004-x64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-U20.yaml@main
+
+  Ubuntu2004-x64-OpenVINO:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-U20-OpenVINO.yaml@main
+
+  Ubuntu2204-x64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-U22.yaml@main
+
+  Ubuntu2404-x64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-U24.yaml@main
+
+  Ubuntu2004-x64-CUDA:
+    if: "${{ contains(github.event.pull_request.labels.*.name, 'category: dnn') }} || ${{ contains(github.event.pull_request.labels.*.name, 'category: dnn (onnx)') }}"
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-U20-Cuda.yaml@main
+
+  Windows10-x64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-W10.yaml@main
+
+  Windows10-x64-UWP:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-W10-UWP.yaml@main
+
+  Windows10-ARM64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-W10-ARM64.yaml@main
+
+# Vulkan configuration disabled as Vulkan backend for DNN does not support int/int64 for now
+# Details: https://github.com/opencv/opencv/issues/25110
+#  Windows10-x64-Vulkan:
+#    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-W10-Vulkan.yaml@main
+
+  macOS-ARM64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-macOS-ARM64.yaml@main
+
+#  macOS-ARM64-Vulkan:
+#    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-macOS-ARM64-Vulkan.yaml@main
+
+  macOS-x64:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-macOS-x86_64.yaml@main
+
+  iOS:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-iOS.yaml@main
+
+  Android:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-Android.yaml@main
+
+  TIM-VX:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-timvx-backend-tests-4.x.yml@main
+
+  docs:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-docs.yaml@main
+
+  Linux-RISC-V-Clang:
+    uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-5.x-RISCV.yaml@main
@@ -21,7 +21,6 @@ jobs:
                crossbuild-essential-arm64 \
                git \
                cmake \
-                libpython-dev:arm64 \
                libpython3-dev:arm64 \
                python-numpy \
                python3-numpy
@@ -34,10 +33,7 @@ jobs:
      run: |
        mkdir build
        cd build
-        cmake -DPYTHON2_INCLUDE_PATH=/usr/include/python2.7/ \
-              -DPYTHON2_LIBRARIES=/usr/lib/aarch64-linux-gnu/libpython2.7.so \
-              -DPYTHON2_NUMPY_INCLUDE_DIRS=/usr/lib/python2.7/dist-packages/numpy/core/include \
-              -DPYTHON3_INCLUDE_PATH=/usr/include/python3.6m/ \
+        cmake -DPYTHON3_INCLUDE_PATH=/usr/include/python3.6m/ \
              -DPYTHON3_LIBRARIES=/usr/lib/aarch64-linux-gnu/libpython3.6m.so \
              -DPYTHON3_NUMPY_INCLUDE_DIRS=/usr/lib/python3/dist-packages/numpy/core/include \
              -DCMAKE_TOOLCHAIN_FILE=../platforms/linux/aarch64-gnu.toolchain.cmake \
@@ -119,7 +119,7 @@ private: \
 #define TEGRA_BINARYOP(type, op, src1, sz1, src2, sz2, dst, sz, w, h) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, h), \
+    parallel_for_(Range(0, h), \
    TegraGenOp_##op##_Invoker<const type, type>(src1, sz1, src2, sz2, dst, sz, w, h), \
    (w * h) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -154,7 +154,7 @@ TegraUnaryOp_Invoker(bitwiseNot, bitwiseNot)
 #define TEGRA_UNARYOP(type, op, src1, sz1, dst, sz, w, h) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, h), \
+    parallel_for_(Range(0, h), \
    TegraGenOp_##op##_Invoker<const type, type>(src1, sz1, dst, sz, w, h), \
    (w * h) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -254,32 +254,32 @@ TegraGenOp_Invoker(cmpLE, cmpGE, 2, 1, 0, RANGE_DATA(ST, src2_data, src2_step),
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
        ((op) == cv::CMP_EQ) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpEQ_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        ((op) == cv::CMP_NE) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpNE_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        ((op) == cv::CMP_GT) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpGT_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        ((op) == cv::CMP_GE) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpGE_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        ((op) == cv::CMP_LT) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpLT_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        ((op) == cv::CMP_LE) ? \
-        parallel_for_(cv::Range(0, h), \
+        parallel_for_(Range(0, h), \
        TegraGenOp_cmpLE_Invoker<const type, CAROTENE_NS::u8>(src1, sz1, src2, sz2, dst, sz, w, h), \
        (w * h) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
@@ -310,7 +310,7 @@ TegraGenOp_Invoker(cmpLE, cmpGE, 2, 1, 0, RANGE_DATA(ST, src2_data, src2_step),
 #define TEGRA_BINARYOPSCALE(type, op, src1, sz1, src2, sz2, dst, sz, w, h, scales) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, h), \
+    parallel_for_(Range(0, h), \
    TegraGenOp_##op##_Invoker<const type, type>(src1, sz1, src2, sz2, dst, sz, w, h, scales), \
    (w * h) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -332,7 +332,7 @@ TegraBinaryOpScale_Invoker(divf, div, 1, scale)
 #define TEGRA_UNARYOPSCALE(type, op, src1, sz1, dst, sz, w, h, scales) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, h), \
+    parallel_for_(Range(0, h), \
    TegraGenOp_##op##_Invoker<const type, type>(src1, sz1, dst, sz, w, h, scales), \
    (w * h) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -928,17 +928,17 @@ TegraRowOp_Invoker(split4, split4, 1, 4, 0, RANGE_DATA(ST, src1_data, 4*sizeof(S
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
        cn == 2 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_split2_Invoker<const type, type>(src, dst[0], dst[1]), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        cn == 3 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_split3_Invoker<const type, type>(src, dst[0], dst[1], dst[2]), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        cn == 4 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_split4_Invoker<const type, type>(src, dst[0], dst[1], dst[2], dst[3]), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
@@ -990,17 +990,17 @@ TegraRowOp_Invoker(combine4, combine4, 4, 1, 0, RANGE_DATA(ST, src1_data, sizeof
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
        cn == 2 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_combine2_Invoker<const type, type>(src[0], src[1], dst), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        cn == 3 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_combine3_Invoker<const type, type>(src[0], src[1], src[2], dst), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
        cn == 4 ? \
-        parallel_for_(cv::Range(0, len), \
+        parallel_for_(Range(0, len), \
        TegraRowOp_combine4_Invoker<const type, type>(src[0], src[1], src[2], src[3], dst), \
        (len) / static_cast<double>(1<<16)), \
        CV_HAL_ERROR_OK : \
@@ -1033,7 +1033,7 @@ TegraRowOp_Invoker(phase, phase, 2, 1, 1, RANGE_DATA(ST, src1_data, sizeof(CAROT
 #define TEGRA_FASTATAN(y, x, dst, len, angleInDegrees) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, len), \
+    parallel_for_(Range(0, len), \
    TegraRowOp_phase_Invoker<const CAROTENE_NS::f32, CAROTENE_NS::f32>(x, y, dst, angleInDegrees ? 1.0f : M_PI/180), \
    (len) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -1049,7 +1049,7 @@ TegraRowOp_Invoker(magnitude, magnitude, 2, 1, 0, RANGE_DATA(ST, src1_data, size
 #define TEGRA_MAGNITUDE(x, y, dst, len) \
 ( \
    CAROTENE_NS::isSupportedConfiguration() ? \
-    parallel_for_(cv::Range(0, len), \
+    parallel_for_(Range(0, len), \
    TegraRowOp_magnitude_Invoker<const CAROTENE_NS::f32, CAROTENE_NS::f32>(x, y, dst), \
    (len) / static_cast<double>(1<<16)), \
    CV_HAL_ERROR_OK \
@@ -1563,17 +1563,17 @@ TegraCvtColor_Invoker(rgbx2bgrx, rgbx2bgrx, src_data + static_cast<size_t>(range
        scn == 3 ? \
            dcn == 3 ? \
                swapBlue ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgb2bgr_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)), \
                    CV_HAL_ERROR_OK : \
                    CV_HAL_ERROR_NOT_IMPLEMENTED : \
            dcn == 4 ? \
                (swapBlue ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgb2bgrx_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgb2rgbx_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) ), \
                CV_HAL_ERROR_OK : \
@@ -1581,16 +1581,16 @@ TegraCvtColor_Invoker(rgbx2bgrx, rgbx2bgrx, src_data + static_cast<size_t>(range
        scn == 4 ? \
            dcn == 3 ? \
                (swapBlue ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgbx2bgr_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgbx2rgb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) ), \
                CV_HAL_ERROR_OK : \
            dcn == 4 ? \
                swapBlue ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgbx2bgrx_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)), \
                    CV_HAL_ERROR_OK : \
@@ -1613,19 +1613,19 @@ TegraCvtColor_Invoker(rgbx2rgb565, rgbx2rgb565, src_data + static_cast<size_t>(r
    greenBits == 6 && CAROTENE_NS::isSupportedConfiguration() ? \
        scn == 3 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgb2bgr565_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgb2rgb565_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
        scn == 4 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgbx2bgr565_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgbx2rgb565_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
@@ -1646,19 +1646,19 @@ TegraCvtColor_Invoker(bgrx2gray, bgrx2gray, CAROTENE_NS::COLOR_SPACE_BT601, src_
    depth == CV_8U && CAROTENE_NS::isSupportedConfiguration() ? \
        scn == 3 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgb2gray_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_bgr2gray_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
        scn == 4 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgbx2gray_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_bgrx2gray_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
@@ -1674,12 +1674,12 @@ TegraCvtColor_Invoker(gray2rgbx, gray2rgbx, src_data + static_cast<size_t>(range
 ( \
    depth == CV_8U && CAROTENE_NS::isSupportedConfiguration() ? \
        dcn == 3 ? \
-            parallel_for_(cv::Range(0, height), \
+            parallel_for_(Range(0, height), \
            TegraCvtColor_gray2rgb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
            (width * height) / static_cast<double>(1<<16)), \
            CV_HAL_ERROR_OK : \
        dcn == 4 ? \
-            parallel_for_(cv::Range(0, height), \
+            parallel_for_(Range(0, height), \
            TegraCvtColor_gray2rgbx_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
            (width * height) / static_cast<double>(1<<16)), \
            CV_HAL_ERROR_OK : \
@@ -1700,19 +1700,19 @@ TegraCvtColor_Invoker(bgrx2ycrcb, bgrx2ycrcb, src_data + static_cast<size_t>(ran
    isCbCr && depth == CV_8U && CAROTENE_NS::isSupportedConfiguration() ? \
        scn == 3 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgb2ycrcb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_bgr2ycrcb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
        scn == 4 ? \
            (swapBlue ? \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_rgbx2ycrcb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) : \
-                parallel_for_(cv::Range(0, height), \
+                parallel_for_(Range(0, height), \
                TegraCvtColor_bgrx2ycrcb_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
@@ -1742,34 +1742,34 @@ TegraCvtColor_Invoker(bgrx2hsvf, bgrx2hsv, src_data + static_cast<size_t>(range.
        scn == 3 ? \
            (swapBlue ? \
                isFullRange ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgb2hsvf_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgb2hsv_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
                isFullRange ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_bgr2hsvf_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_bgr2hsv_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
        scn == 4 ? \
            (swapBlue ? \
                isFullRange ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgbx2hsvf_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_rgbx2hsv_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
                isFullRange ? \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_bgrx2hsvf_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) : \
-                    parallel_for_(cv::Range(0, height), \
+                    parallel_for_(Range(0, height), \
                    TegraCvtColor_bgrx2hsv_Invoker(src_data, src_step, dst_data, dst_step, width, height), \
                    (width * height) / static_cast<double>(1<<16)) ), \
            CV_HAL_ERROR_OK : \
@@ -359,7 +359,7 @@ namespace CAROTENE_NS {

    /*
        For each point `p` within `size`, do:
-        dst[p] = src0[p] * scale / src1[p]
+        dst[p] = src0[p] * scale / src1[p] 

        NOTE: ROUND_TO_ZERO convert policy is used
    */
@@ -420,7 +420,7 @@ namespace CAROTENE_NS {

    /*
        For each point `p` within `size`, do:
-        dst[p] = scale / src[p]
+        dst[p] = scale / src[p] 

        NOTE: ROUND_TO_ZERO convert policy is used
    */
@@ -1040,7 +1040,7 @@ namespace CAROTENE_NS {
                        s32 maxVal, size_t * maxLocPtr, s32 & maxLocCount, s32 maxLocCapacity);

    /*
-        Among each pixel `p` within `src` find min and max values and its first occurrences
+        Among each pixel `p` within `src` find min and max values and its first occurences
    */
    void minMaxLoc(const Size2D &size,
                   const s8 * srcBase, ptrdiff_t srcStride,
@@ -378,7 +378,7 @@ void extract4(const Size2D &size,
                                  vst1q_##sgn##bits(dst1 + d1j, vals.v4.val[3]); \
                              }

-#endif
+#endif 

 #define SPLIT4ALPHA(sgn,bits) void split4(const Size2D &_size,                                          \
                                          const sgn##bits * srcBase, ptrdiff_t srcStride,               \
@@ -0,0 +1,49 @@
+# ----------------------------------------------------------------------------
+#  CMake file for opencv_lapack. See root CMakeLists.txt
+#
+# ----------------------------------------------------------------------------
+project(clapack)
+
+# TODO: extract it from sources somehow
+set(CLAPACK_VERSION "3.9.0" PARENT_SCOPE)
+
+include_directories("${CMAKE_CURRENT_SOURCE_DIR}/include")
+
+# The .cpp files:
+file(GLOB lapack_srcs src/*.c)
+file(GLOB runtime_srcs runtime/*.c)
+file(GLOB lib_hdrs include/*.h)
+
+# ----------------------------------------------------------------------------------
+# 				Define the library target:
+# ----------------------------------------------------------------------------------
+
+set(the_target "libclapack")
+
+add_library(${the_target} STATIC ${lapack_srcs} ${runtime_srcs} ${lib_hdrs})
+
+ocv_warnings_disable(CMAKE_C_FLAGS -Wno-parentheses -Wno-uninitialized -Wno-array-bounds
+    -Wno-implicit-function-declaration -Wno-unused -Wunused-parameter -Wstringop-truncation
+    -Wtautological-negation-compare) # gcc/clang warnings
+ocv_warnings_disable(CMAKE_C_FLAGS /wd4244 /wd4554 /wd4723 /wd4819) # visual studio warnings
+
+set_target_properties(${the_target}
+  PROPERTIES OUTPUT_NAME ${the_target}
+  DEBUG_POSTFIX "${OPENCV_DEBUG_POSTFIX}"
+  COMPILE_PDB_NAME ${the_target}
+  COMPILE_PDB_NAME_DEBUG "${the_target}${OPENCV_DEBUG_POSTFIX}"
+  ARCHIVE_OUTPUT_DIRECTORY ${3P_LIBRARY_OUTPUT_PATH}
+  )
+
+set(CLAPACK_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/include" PARENT_SCOPE)
+set(CLAPACK_LIBRARIES ${the_target} PARENT_SCOPE)
+
+if(ENABLE_SOLUTION_FOLDERS)
+  set_target_properties(${the_target} PROPERTIES FOLDER "3rdparty")
+endif()
+
+if(NOT BUILD_SHARED_LIBS)
+  ocv_install_target(${the_target} EXPORT OpenCVModules ARCHIVE DESTINATION ${OPENCV_3P_LIB_INSTALL_PATH} COMPONENT dev)
+endif()
+
+ocv_install_3rdparty_licenses(clapack lapack_LICENSE)
@@ -0,0 +1,102 @@
+#ifndef __CBLAS_H__
+#define __CBLAS_H__
+
+/* most of the stuff is in lapacke.h */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct lapack_complex
+{
+    float r, i;
+} lapack_complex;
+
+typedef struct lapack_doublecomplex
+{
+    double r, i;
+} lapack_doublecomplex;
+
+typedef enum {CblasRowMajor=101, CblasColMajor=102} CBLAS_LAYOUT;
+typedef enum {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
+
+void cblas_xerbla(const CBLAS_LAYOUT layout, int info,
+                  const char *rout, const char *form, ...);
+
+void cblas_sgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA,
+                 CBLAS_TRANSPOSE TransB, const int M, const int N,
+                 const int K, const float alpha, const float *A,
+                 const int lda, const float *B, const int ldb,
+                 const float beta, float *C, const int ldc);
+
+void cblas_dgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA,
+              CBLAS_TRANSPOSE TransB, const int M, const int N,
+              const int K, const double alpha, const double *A,
+              const int lda, const double *B, const int ldb,
+              const double beta, double *C, const int ldc);
+
+void cblas_cgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA,
+               CBLAS_TRANSPOSE TransB, const int M, const int N,
+               const int K, const void *alpha, const void *A,
+               const int lda, const void *B, const int ldb,
+               const void *beta, void *C, const int ldc);
+
+void cblas_zgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA,
+                CBLAS_TRANSPOSE TransB, const int M, const int N,
+                const int K, const void *alpha, const void *A,
+                const int lda, const void *B, const int ldb,
+                const void *beta, void *C, const int ldc);
+
+int xerbla_(char *, int *);
+int lsame_(char *, char *);
+double slamch_(char* cmach);
+double slamc3_(float *a, float *b);
+double dlamch_(char* cmach);
+double dlamc3_(double *a, double *b);
+
+int dgels_(char *trans, int *m, int *n, int *nrhs, double *a,
+	int *lda, double *b, int *ldb, double *work, int *lwork, int *info);
+
+int dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv,
+	double *b, int *ldb, int *info);
+
+int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv,
+	int *info);
+
+int dposv_(char *uplo, int *n, int *nrhs, double *a, int *
+	lda, double *b, int *ldb, int *info);
+
+int dpotrf_(char *uplo, int *n, double *a, int *lda, int *
+	info);
+
+int sgels_(char *trans, int *m, int *n, int *nrhs, float *a,
+	int *lda, float *b, int *ldb, float *work, int *lwork, int *info);
+
+int sgeev_(char *jobvl, char *jobvr, int *n, float *a, int *
+	lda, float *wr, float *wi, float *vl, int *ldvl, float *vr, int *
+	ldvr, float *work, int *lwork, int *info);
+
+int sgeqrf_(int *m, int *n, float *a, int *lda, float *tau,
+	float *work, int *lwork, int *info);
+
+int sgesv_(int *n, int *nrhs, float *a, int *lda, int *ipiv,
+	float *b, int *ldb, int *info);
+
+int sgetrf_(int *m, int *n, float *a, int *lda, int *ipiv,
+	int *info);
+
+int sposv_(char *uplo, int *n, int *nrhs, float *a, int *
+	lda, float *b, int *ldb, int *info);
+
+int spotrf_(char *uplo, int *n, float *a, int *lda, int *
+	info);
+
+int sgesdd_(char *jobz, int *m, int *n, float *a, int *lda,
+	float *s, float *u, int *ldu, float *vt, int *ldvt, float *work,
+	int *lwork, int *iwork, int *info);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* __CBLAS_H__ */
@@ -0,0 +1,129 @@
+/* f2c.h  --  Standard Fortran to C header file */
+
+/**  barf  [ba:rf]  2.  "He suggested using FORTRAN, and everybody barfed."
+
+	- From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */
+
+#ifndef __F2C_H__
+#define __F2C_H__
+
+#include <assert.h>
+#include <math.h>
+#include <ctype.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "cblas.h"
+#include "lapack.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#undef complex
+
+typedef int integer;
+typedef unsigned int uinteger;
+typedef char *address;
+typedef short int shortint;
+typedef float real;
+typedef double doublereal;
+typedef lapack_complex complex;
+typedef lapack_doublecomplex doublecomplex;
+typedef int logical;
+typedef short int shortlogical;
+typedef char logical1;
+typedef char integer1;
+
+#define TRUE_ (1)
+#define FALSE_ (0)
+
+#ifndef abs
+#define abs(x) ((x) >= 0 ? (x) : -(x))
+#endif
+#define dabs(x) (double)abs(x)
+#ifndef min
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#endif
+#ifndef max
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+#endif
+#define dmin(a,b) (double)min(a,b)
+#define dmax(a,b) (double)max(a,b)
+#define bit_test(a,b)	((a) >> (b) & 1)
+#define bit_clear(a,b)	((a) & ~((uinteger)1 << (b)))
+#define bit_set(a,b)	((a) |  ((uinteger)1 << (b)))
+
+static __inline double r_lg10(float *x)
+{
+    return 0.43429448190325182765*log(*x);
+}
+
+static __inline double d_lg10(double *x)
+{
+    return 0.43429448190325182765*log(*x);
+}
+
+static __inline double d_sign(double *a, double *b)
+{
+    double x = fabs(*a);
+    return *b >= 0 ? x : -x;
+}
+
+static __inline double r_sign(float *a, float *b)
+{
+    double x = fabs((double)*a);
+    return *b >= 0 ? x : -x;
+}
+
+static __inline int i_nint(float *x)
+{
+    return (int)(*x >= 0 ? floor(*x + .5) : -floor(.5 - *x));
+}
+
+int pow_ii(int *ap, int *bp);
+double pow_di(double *ap, int *bp);
+static __inline double pow_ri(float *ap, int *bp)
+{
+    double apd = *ap;
+    return pow_di(&apd, bp);
+}
+static __inline double pow_dd(double *ap, double *bp)
+{
+    return pow(*ap, *bp);
+}
+
+static __inline void d_cnjg(doublecomplex *r, doublecomplex *z)
+{
+	double zi = z->i;
+	r->r = z->r;
+	r->i = -zi;
+}
+
+static __inline void r_cnjg(complex *r, complex *z)
+{
+	float zi = z->i;
+	r->r = z->r;
+	r->i = -zi;
+}
+
+static __inline int s_copy(char *a, char *b, int maxlen)
+{
+    strncpy(a, b, maxlen);
+    a[maxlen] = '\0';
+    return 0;
+}
+
+int s_cat(char *lp, char **rpp, int* rnp, int *np);
+int s_cmp(char *a0, char *b0);
+static __inline int i_len(char* s)
+{
+    return (int)strlen(s);
+}
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
@@ -0,0 +1,386 @@
+// this is auto-generated header for Lapack subset
+#ifndef __CLAPACK_H__
+#define __CLAPACK_H__
+
+#include "cblas.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int cgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, lapack_complex *alpha, lapack_complex *a, int *lda, lapack_complex *b, int *ldb,
+	lapack_complex *beta, lapack_complex *c__, int *ldc);
+
+int daxpy_(int *n, double *da, double *dx, int *incx, double
+	*dy, int *incy);
+
+int dbdsdc_(char *uplo, char *compq, int *n, double *d__,
+	double *e, double *u, int *ldu, double *vt, int *ldvt, double *q, int
+	*iq, double *work, int *iwork, int *info);
+
+int dbdsqr_(char *uplo, int *n, int *ncvt, int *nru, int *
+	ncc, double *d__, double *e, double *vt, int *ldvt, double *u, int *
+	ldu, double *c__, int *ldc, double *work, int *info);
+
+int dcombssq_(double *v1, double *v2);
+
+int dcopy_(int *n, double *dx, int *incx, double *dy, int *
+	incy);
+
+double ddot_(int *n, double *dx, int *incx, double *dy, int *incy);
+
+int dgebak_(char *job, char *side, int *n, int *ilo, int *
+	ihi, double *scale, int *m, double *v, int *ldv, int *info);
+
+int dgebal_(char *job, int *n, double *a, int *lda, int *ilo,
+	 int *ihi, double *scale, int *info);
+
+int dgebd2_(int *m, int *n, double *a, int *lda, double *d__,
+	 double *e, double *tauq, double *taup, double *work, int *info);
+
+int dgebrd_(int *m, int *n, double *a, int *lda, double *d__,
+	 double *e, double *tauq, double *taup, double *work, int *lwork, int
+	*info);
+
+int dgeev_(char *jobvl, char *jobvr, int *n, double *a, int *
+	lda, double *wr, double *wi, double *vl, int *ldvl, double *vr, int *
+	ldvr, double *work, int *lwork, int *info);
+
+int dgehd2_(int *n, int *ilo, int *ihi, double *a, int *lda,
+	double *tau, double *work, int *info);
+
+int dgehrd_(int *n, int *ilo, int *ihi, double *a, int *lda,
+	double *tau, double *work, int *lwork, int *info);
+
+int dgelq2_(int *m, int *n, double *a, int *lda, double *tau,
+	 double *work, int *info);
+
+int dgelqf_(int *m, int *n, double *a, int *lda, double *tau,
+	 double *work, int *lwork, int *info);
+
+int dgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, double *alpha, double *a, int *lda, double *b, int *ldb, double *
+	beta, double *c__, int *ldc);
+
+int dgemv_(char *trans, int *m, int *n, double *alpha,
+	double *a, int *lda, double *x, int *incx, double *beta, double *y,
+	int *incy);
+
+int dgeqr2_(int *m, int *n, double *a, int *lda, double *tau,
+	 double *work, int *info);
+
+int dgeqrf_(int *m, int *n, double *a, int *lda, double *tau,
+	 double *work, int *lwork, int *info);
+
+int dger_(int *m, int *n, double *alpha, double *x, int *
+	incx, double *y, int *incy, double *a, int *lda);
+
+int dgesdd_(char *jobz, int *m, int *n, double *a, int *lda,
+	double *s, double *u, int *ldu, double *vt, int *ldvt, double *work,
+	int *lwork, int *iwork, int *info);
+
+int dhseqr_(char *job, char *compz, int *n, int *ilo, int *
+	ihi, double *h__, int *ldh, double *wr, double *wi, double *z__, int *
+	ldz, double *work, int *lwork, int *info);
+
+int disnan_(double *din);
+
+// "small" is a macro defined in Windows headers: https://stackoverflow.com/a/27794577
+#ifdef small
+#undef small
+#endif
+
+int dlabad_(double *small, double *large);
+
+int dlabrd_(int *m, int *n, int *nb, double *a, int *lda,
+	double *d__, double *e, double *tauq, double *taup, double *x, int *
+	ldx, double *y, int *ldy);
+
+int dlacpy_(char *uplo, int *m, int *n, double *a, int *lda,
+	double *b, int *ldb);
+
+int dladiv1_(double *a, double *b, double *c__, double *d__,
+	double *p, double *q);
+
+double dladiv2_(double *a, double *b, double *c__, double *d__, double *r__,
+	double *t);
+
+int dladiv_(double *a, double *b, double *c__, double *d__,
+	double *p, double *q);
+
+int dlaed6_(int *kniter, int *orgati, double *rho, double *
+	d__, double *z__, double *finit, double *tau, int *info);
+
+int dlaexc_(int *wantq, int *n, double *t, int *ldt, double *
+	q, int *ldq, int *j1, int *n1, int *n2, double *work, int *info);
+
+int dlahqr_(int *wantt, int *wantz, int *n, int *ilo, int *
+	ihi, double *h__, int *ldh, double *wr, double *wi, int *iloz, int *
+	ihiz, double *z__, int *ldz, int *info);
+
+int dlahr2_(int *n, int *k, int *nb, double *a, int *lda,
+	double *tau, double *t, int *ldt, double *y, int *ldy);
+
+int dlaisnan_(double *din1, double *din2);
+
+int dlaln2_(int *ltrans, int *na, int *nw, double *smin,
+	double *ca, double *a, int *lda, double *d1, double *d2, double *b,
+	int *ldb, double *wr, double *wi, double *x, int *ldx, double *scale,
+	double *xnorm, int *info);
+
+int dlamrg_(int *n1, int *n2, double *a, int *dtrd1, int *
+	dtrd2, int *index);
+
+double dlange_(char *norm, int *m, int *n, double *a, int *lda, double *work);
+
+double dlanst_(char *norm, int *n, double *d__, double *e);
+
+int dlanv2_(double *a, double *b, double *c__, double *d__,
+	double *rt1r, double *rt1i, double *rt2r, double *rt2i, double *cs,
+	double *sn);
+
+double dlapy2_(double *x, double *y);
+
+int dlaqr0_(int *wantt, int *wantz, int *n, int *ilo, int *
+	ihi, double *h__, int *ldh, double *wr, double *wi, int *iloz, int *
+	ihiz, double *z__, int *ldz, double *work, int *lwork, int *info);
+
+int dlaqr1_(int *n, double *h__, int *ldh, double *sr1,
+	double *si1, double *sr2, double *si2, double *v);
+
+int dlaqr2_(int *wantt, int *wantz, int *n, int *ktop, int *
+	kbot, int *nw, double *h__, int *ldh, int *iloz, int *ihiz, double *
+	z__, int *ldz, int *ns, int *nd, double *sr, double *si, double *v,
+	int *ldv, int *nh, double *t, int *ldt, int *nv, double *wv, int *
+	ldwv, double *work, int *lwork);
+
+int dlaqr3_(int *wantt, int *wantz, int *n, int *ktop, int *
+	kbot, int *nw, double *h__, int *ldh, int *iloz, int *ihiz, double *
+	z__, int *ldz, int *ns, int *nd, double *sr, double *si, double *v,
+	int *ldv, int *nh, double *t, int *ldt, int *nv, double *wv, int *
+	ldwv, double *work, int *lwork);
+
+int dlaqr4_(int *wantt, int *wantz, int *n, int *ilo, int *
+	ihi, double *h__, int *ldh, double *wr, double *wi, int *iloz, int *
+	ihiz, double *z__, int *ldz, double *work, int *lwork, int *info);
+
+int dlaqr5_(int *wantt, int *wantz, int *kacc22, int *n, int
+	*ktop, int *kbot, int *nshfts, double *sr, double *si, double *h__,
+	int *ldh, int *iloz, int *ihiz, double *z__, int *ldz, double *v, int
+	*ldv, double *u, int *ldu, int *nv, double *wv, int *ldwv, int *nh,
+	double *wh, int *ldwh);
+
+int dlarf_(char *side, int *m, int *n, double *v, int *incv,
+	double *tau, double *c__, int *ldc, double *work);
+
+int dlarfb_(char *side, char *trans, char *direct, char *
+	storev, int *m, int *n, int *k, double *v, int *ldv, double *t, int *
+	ldt, double *c__, int *ldc, double *work, int *ldwork);
+
+int dlarfg_(int *n, double *alpha, double *x, int *incx,
+	double *tau);
+
+int dlarft_(char *direct, char *storev, int *n, int *k,
+	double *v, int *ldv, double *tau, double *t, int *ldt);
+
+int dlarfx_(char *side, int *m, int *n, double *v, double *
+	tau, double *c__, int *ldc, double *work);
+
+int dlartg_(double *f, double *g, double *cs, double *sn,
+	double *r__);
+
+int dlas2_(double *f, double *g, double *h__, double *ssmin,
+	double *ssmax);
+
+int dlascl_(char *type__, int *kl, int *ku, double *cfrom,
+	double *cto, int *m, int *n, double *a, int *lda, int *info);
+
+int dlasd0_(int *n, int *sqre, double *d__, double *e,
+	double *u, int *ldu, double *vt, int *ldvt, int *smlsiz, int *iwork,
+	double *work, int *info);
+
+int dlasd1_(int *nl, int *nr, int *sqre, double *d__, double
+	*alpha, double *beta, double *u, int *ldu, double *vt, int *ldvt, int
+	*idxq, int *iwork, double *work, int *info);
+
+int dlasd2_(int *nl, int *nr, int *sqre, int *k, double *d__,
+	 double *z__, double *alpha, double *beta, double *u, int *ldu,
+	double *vt, int *ldvt, double *dsigma, double *u2, int *ldu2, double *
+	vt2, int *ldvt2, int *idxp, int *idx, int *idxc, int *idxq, int *
+	coltyp, int *info);
+
+int dlasd3_(int *nl, int *nr, int *sqre, int *k, double *d__,
+	 double *q, int *ldq, double *dsigma, double *u, int *ldu, double *u2,
+	 int *ldu2, double *vt, int *ldvt, double *vt2, int *ldvt2, int *idxc,
+	 int *ctot, double *z__, int *info);
+
+int dlasd4_(int *n, int *i__, double *d__, double *z__,
+	double *delta, double *rho, double *sigma, double *work, int *info);
+
+int dlasd5_(int *i__, double *d__, double *z__, double *
+	delta, double *rho, double *dsigma, double *work);
+
+int dlasd6_(int *icompq, int *nl, int *nr, int *sqre, double
+	*d__, double *vf, double *vl, double *alpha, double *beta, int *idxq,
+	int *perm, int *givptr, int *givcol, int *ldgcol, double *givnum, int
+	*ldgnum, double *poles, double *difl, double *difr, double *z__, int *
+	k, double *c__, double *s, double *work, int *iwork, int *info);
+
+int dlasd7_(int *icompq, int *nl, int *nr, int *sqre, int *k,
+	 double *d__, double *z__, double *zw, double *vf, double *vfw,
+	double *vl, double *vlw, double *alpha, double *beta, double *dsigma,
+	int *idx, int *idxp, int *idxq, int *perm, int *givptr, int *givcol,
+	int *ldgcol, double *givnum, int *ldgnum, double *c__, double *s, int
+	*info);
+
+int dlasd8_(int *icompq, int *k, double *d__, double *z__,
+	double *vf, double *vl, double *difl, double *difr, int *lddifr,
+	double *dsigma, double *work, int *info);
+
+int dlasda_(int *icompq, int *smlsiz, int *n, int *sqre,
+	double *d__, double *e, double *u, int *ldu, double *vt, int *k,
+	double *difl, double *difr, double *z__, double *poles, int *givptr,
+	int *givcol, int *ldgcol, int *perm, double *givnum, double *c__,
+	double *s, double *work, int *iwork, int *info);
+
+int dlasdq_(char *uplo, int *sqre, int *n, int *ncvt, int *
+	nru, int *ncc, double *d__, double *e, double *vt, int *ldvt, double *
+	u, int *ldu, double *c__, int *ldc, double *work, int *info);
+
+int dlasdt_(int *n, int *lvl, int *nd, int *inode, int *
+	ndiml, int *ndimr, int *msub);
+
+int dlaset_(char *uplo, int *m, int *n, double *alpha,
+	double *beta, double *a, int *lda);
+
+int dlasq1_(int *n, double *d__, double *e, double *work,
+	int *info);
+
+int dlasq2_(int *n, double *z__, int *info);
+
+int dlasq3_(int *i0, int *n0, double *z__, int *pp, double *
+	dmin__, double *sigma, double *desig, double *qmax, int *nfail, int *
+	iter, int *ndiv, int *ieee, int *ttype, double *dmin1, double *dmin2,
+	double *dn, double *dn1, double *dn2, double *g, double *tau);
+
+int dlasq4_(int *i0, int *n0, double *z__, int *pp, int *
+	n0in, double *dmin__, double *dmin1, double *dmin2, double *dn,
+	double *dn1, double *dn2, double *tau, int *ttype, double *g);
+
+int dlasq5_(int *i0, int *n0, double *z__, int *pp, double *
+	tau, double *sigma, double *dmin__, double *dmin1, double *dmin2,
+	double *dn, double *dnm1, double *dnm2, int *ieee, double *eps);
+
+int dlasq6_(int *i0, int *n0, double *z__, int *pp, double *
+	dmin__, double *dmin1, double *dmin2, double *dn, double *dnm1,
+	double *dnm2);
+
+int dlasr_(char *side, char *pivot, char *direct, int *m,
+	int *n, double *c__, double *s, double *a, int *lda);
+
+int dlasrt_(char *id, int *n, double *d__, int *info);
+
+int dlassq_(int *n, double *x, int *incx, double *scale,
+	double *sumsq);
+
+int dlasv2_(double *f, double *g, double *h__, double *ssmin,
+	 double *ssmax, double *snr, double *csr, double *snl, double *csl);
+
+int dlasy2_(int *ltranl, int *ltranr, int *isgn, int *n1,
+	int *n2, double *tl, int *ldtl, double *tr, int *ldtr, double *b, int
+	*ldb, double *scale, double *x, int *ldx, double *xnorm, int *info);
+
+double dnrm2_(int *n, double *x, int *incx);
+
+int dorg2r_(int *m, int *n, int *k, double *a, int *lda,
+	double *tau, double *work, int *info);
+
+int dorgbr_(char *vect, int *m, int *n, int *k, double *a,
+	int *lda, double *tau, double *work, int *lwork, int *info);
+
+int dorghr_(int *n, int *ilo, int *ihi, double *a, int *lda,
+	double *tau, double *work, int *lwork, int *info);
+
+int dorgl2_(int *m, int *n, int *k, double *a, int *lda,
+	double *tau, double *work, int *info);
+
+int dorglq_(int *m, int *n, int *k, double *a, int *lda,
+	double *tau, double *work, int *lwork, int *info);
+
+int dorgqr_(int *m, int *n, int *k, double *a, int *lda,
+	double *tau, double *work, int *lwork, int *info);
+
+int dorm2r_(char *side, char *trans, int *m, int *n, int *k,
+	double *a, int *lda, double *tau, double *c__, int *ldc, double *work,
+	 int *info);
+
+int dormbr_(char *vect, char *side, char *trans, int *m, int
+	*n, int *k, double *a, int *lda, double *tau, double *c__, int *ldc,
+	double *work, int *lwork, int *info);
+
+int dormhr_(char *side, char *trans, int *m, int *n, int *
+	ilo, int *ihi, double *a, int *lda, double *tau, double *c__, int *
+	ldc, double *work, int *lwork, int *info);
+
+int dorml2_(char *side, char *trans, int *m, int *n, int *k,
+	double *a, int *lda, double *tau, double *c__, int *ldc, double *work,
+	 int *info);
+
+int dormlq_(char *side, char *trans, int *m, int *n, int *k,
+	double *a, int *lda, double *tau, double *c__, int *ldc, double *work,
+	 int *lwork, int *info);
+
+int dormqr_(char *side, char *trans, int *m, int *n, int *k,
+	double *a, int *lda, double *tau, double *c__, int *ldc, double *work,
+	 int *lwork, int *info);
+
+int drot_(int *n, double *dx, int *incx, double *dy, int *
+	incy, double *c__, double *s);
+
+int dscal_(int *n, double *da, double *dx, int *incx);
+
+int dswap_(int *n, double *dx, int *incx, double *dy, int *
+	incy);
+
+int dtrevc3_(char *side, char *howmny, int *select, int *n,
+	double *t, int *ldt, double *vl, int *ldvl, double *vr, int *ldvr,
+	int *mm, int *m, double *work, int *lwork, int *info);
+
+int dtrexc_(char *compq, int *n, double *t, int *ldt, double
+	*q, int *ldq, int *ifst, int *ilst, double *work, int *info);
+
+int dtrmm_(char *side, char *uplo, char *transa, char *diag,
+	int *m, int *n, double *alpha, double *a, int *lda, double *b, int *
+	ldb);
+
+int dtrmv_(char *uplo, char *trans, char *diag, int *n,
+	double *a, int *lda, double *x, int *incx);
+
+int idamax_(int *n, double *dx, int *incx);
+
+int ieeeck_(int *ispec, float *zero, float *one);
+
+int iladlc_(int *m, int *n, double *a, int *lda);
+
+int iladlr_(int *m, int *n, double *a, int *lda);
+
+int ilaenv_(int *ispec, char *name__, char *opts, int *n1, int *n2, int *n3,
+	int *n4);
+
+int iparmq_(int *ispec, char *name__, char *opts, int *n, int *ilo, int *ihi,
+	int *lwork);
+
+int sgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, float *alpha, float *a, int *lda, float *b, int *ldb, float *beta,
+	float *c__, int *ldc);
+
+int zgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, lapack_doublecomplex *alpha, lapack_doublecomplex *a, int *lda, lapack_doublecomplex *b,
+	 int *ldb, lapack_doublecomplex *beta, lapack_doublecomplex *c__, int *ldc);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
@@ -0,0 +1,48 @@
+Copyright (c) 1992-2017 The University of Tennessee and The University
+                        of Tennessee Research Foundation.  All rights
+                        reserved.
+Copyright (c) 2000-2017 The University of California Berkeley. All
+                        rights reserved.
+Copyright (c) 2006-2017 The University of Colorado Denver.  All rights
+                        reserved.
+
+$COPYRIGHT$
+
+Additional copyrights may follow
+
+$HEADER$
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer listed
+  in this license in the documentation and/or other materials
+  provided with the distribution.
+
+- Neither the name of the copyright holders nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+The copyright holders provide no reassurances that the source code
+provided does not infringe any patent, copyright, or any other
+intellectual property rights of third parties.  The copyright holders
+disclaim any liability to any recipient for claims brought against
+recipient by any third party for infringement of that parties
+intellectual property rights.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@@ -0,0 +1,272 @@
+appdoc = """
+    This is generator of CLapack subset.
+    The usage:
+
+    1. Make sure you have the special version of f2c installed.
+       Grab it from https://github.com/vpisarev/f2c/tree/for_lapack.
+    2. Download fresh version of Lapack from
+       https://github.com/Reference-LAPACK/lapack.
+       You may choose some specific version or the latest snapshot.
+    3. If necessary, edit "roots" and "banlist" variables in this script, specify the needed and unneeded functions
+    4. From within a working directory run
+
+       $ python3 <opencv_root>/3rdparty/clapack/make_clapack.py <lapack_root>
+       or
+       $ F2C=<path_to_custom_f2c> python3 <opencv_root>/3rdparty/clapack/make_clapack.py <lapack_root>
+
+       it will generate "new_clapack" directory with "include" and "src" subdirectories.
+    5. erase opencv/3rdparty/clapack/src and replace it with new_clapack/src.
+    6. copy new_clapack/include/lapack.h to opencv/3rdparty/clapack/include.
+    7. optionally, edit opencv/3rdparty/clapack/CMakeLists.txt and update CLAPACK_VERSION as needed.
+
+    This is it. Now build it and enjoy.
+"""
+
+import glob, re, os, shutil, subprocess, sys
+
+roots = ["cgemm_", "dgemm_", "sgemm_", "zgemm_",
+         "dgeev_", "dgesdd_", 
+         #"dsyevr_",
+         #"dgesv_", "dgetrf_", "dposv_", "dpotrf_", "dgels_", "dgeqrf_",
+         #"sgesv_", "sgetrf_", "sposv_", "spotrf_", "sgels_", "sgeqrf_"
+         ]
+banlist = ["slamch_", "slamc3_", "dlamch_", "dlamc3_", "lsame_", "xerbla_"]
+
+if len(sys.argv) < 2:
+    print(appdoc)
+    sys.exit(0)
+
+lapack_root = sys.argv[1]
+dst_path = "."
+
+def error(msg):
+    print ("error: " + msg)
+    sys.exit(0)
+
+def file2fun(fname):
+    return (os.path.basename(fname)[:-2]).upper()
+
+def print_graph(m):
+    for (k, neighbors) in sorted(m.items()):
+        print (k + " : " + ", ".join(sorted(list(neighbors))))
+
+blas_path = os.path.join(lapack_root, "BLAS/SRC")
+lapack_path = os.path.join(lapack_root, "SRC")
+
+roots = [f[:-1].upper() for f in roots]
+banlist = [f[:-1].upper() for f in banlist]
+
+def fun2file(func):
+    filename = func.lower() + ".f"
+    blas_loc = blas_path + "/" + filename
+    lapack_loc = lapack_path + "/" + filename
+    if os.path.exists(blas_loc):
+        return blas_loc
+    elif os.path.exists(lapack_loc):
+        return lapack_loc
+    else:
+        error("neither %s nor %s exist" % (blas_loc, lapack_loc))
+
+all_files = glob.glob(blas_path + "/*.f") + glob.glob(lapack_path + "/*.f")
+all_funcs = [file2fun(fname) for fname in all_files]
+all_funcs_set = set(all_funcs).difference(set(banlist))
+all_funcs = sorted(list(all_funcs_set))
+
+func_deps = {}
+
+#print all_funcs
+
+words_regexp = re.compile(r'\w+')
+
+def scan_deps(func):
+    global func_deps
+    if func in func_deps:
+        return
+    func_deps[func] = set([]) # to avoid possibly infinite recursion
+    f = open(fun2file(func), 'rt')
+    deps = []
+    external_mode = False
+    for l in f.readlines():
+        if l.startswith('*'):
+            continue
+        l = l.strip().upper()
+        if l.startswith('EXTERNAL '):
+            external_mode = True
+        elif l.startswith('$') and external_mode:
+            pass
+        else:
+            external_mode = False
+        if not external_mode:
+            continue
+        for w in words_regexp.findall(l):
+            if w in all_funcs_set:
+                deps.append(w)
+    f.close()
+    # remove func from its dependencies
+    deps = set(deps).difference(set([func]))
+    func_deps[func] = deps
+    for d in deps:
+        scan_deps(d)
+
+for r in roots:
+    scan_deps(r)
+
+selected_funcs = sorted(func_deps.keys())
+print ("total files before amalgamation: %d" % len(selected_funcs))
+
+inv_deps = {}
+for func in selected_funcs:
+    inv_deps[func] = set([])
+
+for (func, deps) in func_deps.items():
+    for d in deps:
+        inv_deps[d] = inv_deps[d].union(set([func]))
+
+#print_graph(inv_deps)
+
+func_home = {}
+for func in selected_funcs:
+    func_home[func] = func
+
+def get_home0(func, func0):
+    used_by = inv_deps[func]
+    if len(used_by) == 1:
+        p = list(used_by)[0]
+        if p != func and p != func0:
+            return get_home0(p, func0)
+        return func
+    return func
+
+# try to merge some files
+for func in selected_funcs:
+    func_home[func] = get_home0(func, func)
+
+# try to merge some files even more
+for iters in range(100):
+    homes_changed = False
+    for (func, used_by) in inv_deps.items():
+        p0 = func_home[func]
+        n = len(used_by)
+        if n == 1:
+            p = list(used_by)[0]
+            p1 = func_home[p]
+            if p1 != p0:
+                func_home[func] = p1
+                homes_changed = True
+            continue
+        elif n > 1:
+            phomes = set([])
+            for p in used_by:
+                phomes.add(func_home[p])
+            if len(phomes) == 1:
+                p1 = list(phomes)[0]
+                if p1 != p0:
+                    func_home[func] = p1
+                    homes_changed = True
+    if not homes_changed:
+        break
+
+res_files = {}
+for (func, h) in func_home.items():
+    elems = res_files.get(h, set([]))
+    elems.add(func)
+    res_files[h] = elems
+
+print ("total files after amalgamation: %d" % len(res_files))
+#print_graph(res_files)
+
+outdir = os.path.join(dst_path, "new_clapack")
+outdir_src = os.path.join(outdir, "src")
+outdir_inc = os.path.join(outdir, "include")
+
+shutil.rmtree(outdir, ignore_errors=True)
+try:
+    os.makedirs(outdir_src)
+except os.error:
+    pass
+try:
+    os.makedirs(outdir_inc)
+except os.error:
+    pass
+
+f2c_appname = os.getenv("F2C", default="f2c")
+print ("f2c used: %s" % f2c_appname)
+
+f2c_getver_cmd = f2c_appname + " -v"
+
+verstr = subprocess.check_output(f2c_getver_cmd.split(' ')).decode("utf-8")
+if "for_lapack" not in verstr:
+    error("invalid version of f2c\n" + appdoc)
+
+f2c_flags = "-ctypes -localconst -no-proto"
+f2c_cmd0 = f2c_appname + " " + f2c_flags
+f2c_cmd1 = f2c_appname + " -hdr none " + f2c_flags
+
+lapack_protos = {}
+extract_fn_regexp = re.compile(r'.+?(\w+)\s*\(')
+
+def extract_proto(func, csrc):
+    global lapack_protos
+    cname = func.lower() + "_"
+    cfname = func.lower() + ".c"
+    regexp_str = r'\n(?:/\* Subroutine \*/\s*)?\w+\s+\w+\s*\((?:.|\n)+?\)[\s\n]*\{'
+    proto_regexp = re.compile(regexp_str)
+    ps = proto_regexp.findall(csrc)
+    for p in ps:
+        n = p.find("*/")
+        if n < 0:
+            n = 0
+        else:
+            n += 2
+        p = p[n:-1].strip() + ";"
+        fns = extract_fn_regexp.findall(p)
+        if len(fns) != 1:
+            error("prototype of function (%s) when analyzing %s cannot be parsed" % (p, cfname))
+        fn = fns[0]
+        if fn not in lapack_protos:
+            p = re.sub(r'\bcomplex\b', 'lapack_complex', p)
+            p = re.sub(r'\bdoublecomplex\b', 'lapack_doublecomplex', p)
+            lapack_protos[fn] = p
+
+for (filename, funcs) in sorted(res_files.items()):
+    out = ""
+    f2c_cmd = f2c_cmd0
+    for func in sorted(list(funcs)):
+        ffilename = fun2file(func)
+        print ("running " + f2c_cmd + " on " + ffilename +  " ...")
+        ffile = open(ffilename, 'rt')
+        delta_out = subprocess.check_output(f2c_cmd.split(' '), stdin=ffile).decode("utf-8")
+        # remove trailing whitespaces
+        delta_out = '\n'.join([l.rstrip() for l in delta_out.split('\n')])
+        extract_proto(func, delta_out)
+        out += delta_out
+        ffile.close()
+        f2c_cmd = f2c_cmd1
+    outname = os.path.join(outdir_src, filename.lower() + ".c")
+    outfile = open(outname, 'wt')
+    outfile.write(out)
+    outfile.close()
+
+proto_hdr = """// this is auto-generated header for Lapack subset
+#ifndef __CLAPACK_H__
+#define __CLAPACK_H__
+
+#include "cblas.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+%s
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+""" % "\n\n".join([p for (n, p) in sorted(lapack_protos.items())])
+
+proto_hdr_fname = os.path.join(outdir_inc, "lapack.h")
+f = open(proto_hdr_fname, 'wt')
+f.write(proto_hdr)
+f.close()
@@ -0,0 +1,289 @@
+#include "f2c.h"
+#include <stdarg.h>
+
+void cblas_cgemm(const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA,
+                 const CBLAS_TRANSPOSE TransB, const int M, const int N,
+                 const int K, const void *alpha, const void  *A,
+                 const int lda, const void  *B, const int ldb,
+                 const void *beta, void  *C, const int ldc)
+{
+   char TA, TB;
+
+   if( layout == CblasColMajor )
+   {
+      if(TransA == CblasTrans) TA='T';
+      else if ( TransA == CblasConjTrans ) TA='C';
+      else if ( TransA == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_cgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+
+      if(TransB == CblasTrans) TB='T';
+      else if ( TransB == CblasConjTrans ) TB='C';
+      else if ( TransB == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 3, "cblas_cgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      cgemm_(&TA, &TB, (int*)&M, (int*)&N, (int*)&K, (complex*)alpha, (complex*)A, (int*)&lda,
+             (complex*)B, (int*)&ldb, (complex*)beta, (complex*)C, (int*)&ldc);
+   }
+   else if (layout == CblasRowMajor)
+   {
+      if(TransA == CblasTrans) TB='T';
+      else if ( TransA == CblasConjTrans ) TB='C';
+      else if ( TransA == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_cgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+      if(TransB == CblasTrans) TA='T';
+      else if ( TransB == CblasConjTrans ) TA='C';
+      else if ( TransB == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_cgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      cgemm_(&TA, &TB, (int*)&N, (int*)&M, (int*)&K, (complex*)alpha, (complex*)B, (int*)&ldb,
+             (complex*)A, (int*)&lda, (complex*)beta, (complex*)C, (int*)&ldc);
+   }
+   else cblas_xerbla(layout, 1, "cblas_cgemm", "Illegal layout setting, %d\n", layout);
+}
+
+void cblas_dgemm(const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA,
+                 const CBLAS_TRANSPOSE TransB, const int M, const int N,
+                 const int K, const double alpha, const double *A,
+                 const int lda, const double  *B, const int ldb,
+                 const double beta, double  *C, const int ldc)
+{
+   char TA, TB;
+
+   if( layout == CblasColMajor )
+   {
+      if(TransA == CblasTrans) TA='T';
+      else if ( TransA == CblasConjTrans ) TA='C';
+      else if ( TransA == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_dgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+
+      if(TransB == CblasTrans) TB='T';
+      else if ( TransB == CblasConjTrans ) TB='C';
+      else if ( TransB == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 3, "cblas_dgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      dgemm_(&TA, &TB, (int*)&M, (int*)&N, (int*)&K, (double*)&alpha, (double*)A, (int*)&lda,
+             (double*)B, (int*)&ldb, (double*)&beta, (double*)C, (int*)&ldc);
+   }
+   else if (layout == CblasRowMajor)
+   {
+      if(TransA == CblasTrans) TB='T';
+      else if ( TransA == CblasConjTrans ) TB='C';
+      else if ( TransA == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_dgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+      if(TransB == CblasTrans) TA='T';
+      else if ( TransB == CblasConjTrans ) TA='C';
+      else if ( TransB == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_dgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      dgemm_(&TA, &TB, (int*)&N, (int*)&M, (int*)&K, (double*)&alpha, (double*)B, (int*)&ldb,
+             (double*)A, (int*)&lda, (double*)&beta, (double*)C, (int*)&ldc);
+   }
+   else cblas_xerbla(layout, 1, "cblas_dgemm", "Illegal layout setting, %d\n", layout);
+}
+
+
+void cblas_sgemm(const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA,
+                 const CBLAS_TRANSPOSE TransB, const int M, const int N,
+                 const int K, const float alpha, const float  *A,
+                 const int lda, const float  *B, const int ldb,
+                 const float beta, float  *C, const int ldc)
+{
+   char TA, TB;
+
+   if( layout == CblasColMajor )
+   {
+      if(TransA == CblasTrans) TA='T';
+      else if ( TransA == CblasConjTrans ) TA='C';
+      else if ( TransA == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_sgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+
+      if(TransB == CblasTrans) TB='T';
+      else if ( TransB == CblasConjTrans ) TB='C';
+      else if ( TransB == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 3, "cblas_sgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      sgemm_(&TA, &TB, (int*)&M, (int*)&N, (int*)&K, (float*)&alpha, (float*)A, (int*)&lda,
+             (float*)B, (int*)&ldb, (float*)&beta, (float*)C, (int*)&ldc);
+   }
+   else if (layout == CblasRowMajor)
+   {
+      if(TransA == CblasTrans) TB='T';
+      else if ( TransA == CblasConjTrans ) TB='C';
+      else if ( TransA == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_sgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+      if(TransB == CblasTrans) TA='T';
+      else if ( TransB == CblasConjTrans ) TA='C';
+      else if ( TransB == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_sgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      sgemm_(&TA, &TB, (int*)&N, (int*)&M, (int*)&K, (float*)&alpha, (float*)B, (int*)&ldb,
+             (float*)A, (int*)&lda, (float*)&beta, (float*)C, (int*)&ldc);
+   }
+   else cblas_xerbla(layout, 1, "cblas_sgemm", "Illegal layout setting, %d\n", layout);
+}
+
+void cblas_zgemm(const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA,
+                 const CBLAS_TRANSPOSE TransB, const int M, const int N,
+                 const int K, const void *alpha, const void  *A,
+                 const int lda, const void  *B, const int ldb,
+                 const void *beta, void  *C, const int ldc)
+{
+   char TA, TB;
+
+   if( layout == CblasColMajor )
+   {
+      if(TransA == CblasTrans) TA='T';
+      else if ( TransA == CblasConjTrans ) TA='C';
+      else if ( TransA == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_zgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+
+      if(TransB == CblasTrans) TB='T';
+      else if ( TransB == CblasConjTrans ) TB='C';
+      else if ( TransB == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 3, "cblas_zgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      zgemm_(&TA, &TB, (int*)&M, (int*)&N, (int*)&K, (doublecomplex*)alpha, (doublecomplex*)A, (int*)&lda,
+             (doublecomplex*)B, (int*)&ldb, (doublecomplex*)beta, (doublecomplex*)C, (int*)&ldc);
+   }
+   else if (layout == CblasRowMajor)
+   {
+      if(TransA == CblasTrans) TB='T';
+      else if ( TransA == CblasConjTrans ) TB='C';
+      else if ( TransA == CblasNoTrans )   TB='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_zgemm", "Illegal TransA setting, %d\n", TransA);
+         return;
+      }
+      if(TransB == CblasTrans) TA='T';
+      else if ( TransB == CblasConjTrans ) TA='C';
+      else if ( TransB == CblasNoTrans )   TA='N';
+      else
+      {
+         cblas_xerbla(layout, 2, "cblas_zgemm", "Illegal TransB setting, %d\n", TransB);
+         return;
+      }
+
+      zgemm_(&TA, &TB, (int*)&N, (int*)&M, (int*)&K, (doublecomplex*)alpha, (doublecomplex*)B, (int*)&ldb,
+             (doublecomplex*)A, (int*)&lda, (doublecomplex*)beta, (doublecomplex*)C, (int*)&ldc);
+   }
+   else cblas_xerbla(layout, 1, "cblas_zgemm", "Illegal layout setting, %d\n", layout);
+}
+
+void cblas_xerbla(const CBLAS_LAYOUT layout, int info, const char *rout, const char *form, ...)
+{
+   extern int RowMajorStrg;
+   char empty[1] = "";
+   va_list argptr;
+
+   va_start(argptr, form);
+
+   if (layout == CblasRowMajor)
+   {
+      if (strstr(rout,"gemm") != 0)
+      {
+         if      (info == 5 ) info =  4;
+         else if (info == 4 ) info =  5;
+         else if (info == 11) info =  9;
+         else if (info == 9 ) info = 11;
+      }
+      else if (strstr(rout,"symm") != 0 || strstr(rout,"hemm") != 0)
+      {
+         if      (info == 5 ) info =  4;
+         else if (info == 4 ) info =  5;
+      }
+      else if (strstr(rout,"trmm") != 0 || strstr(rout,"trsm") != 0)
+      {
+         if      (info == 7 ) info =  6;
+         else if (info == 6 ) info =  7;
+      }
+      else if (strstr(rout,"gemv") != 0)
+      {
+         if      (info == 4)  info = 3;
+         else if (info == 3)  info = 4;
+      }
+      else if (strstr(rout,"gbmv") != 0)
+      {
+         if      (info == 4)  info = 3;
+         else if (info == 3)  info = 4;
+         else if (info == 6)  info = 5;
+         else if (info == 5)  info = 6;
+      }
+      else if (strstr(rout,"ger") != 0)
+      {
+         if      (info == 3) info = 2;
+         else if (info == 2) info = 3;
+         else if (info == 8) info = 6;
+         else if (info == 6) info = 8;
+      }
+      else if ( (strstr(rout,"her2") != 0 || strstr(rout,"hpr2") != 0)
+                 && strstr(rout,"her2k") == 0 )
+      {
+         if      (info == 8) info = 6;
+         else if (info == 6) info = 8;
+      }
+   }
+   if (info)
+      fprintf(stderr, "Parameter %d to routine %s was incorrect\n", info, rout);
+   vfprintf(stderr, form, argptr);
+   va_end(argptr);
+   if (info && !info)
+      xerbla_(empty, &info); /* Force link of our F77 error handler */
+   exit(-1);
+}
@@ -0,0 +1,72 @@
+#include "f2c.h"
+#include <float.h>
+#include <stdio.h>
+
+/* *********************************************************************** */
+
+double dlamc3_(double *a, double *b)
+{
+    /*  -- LAPACK auxiliary routine (version 3.1) -- */
+    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+    /*     November 2006 */
+
+    /*     .. Scalar Arguments .. */
+    /*     .. */
+
+    /*  Purpose */
+    /*  ======= */
+
+    /*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing */
+    /*  the addition of  A  and  B ,  for use in situations where optimizers */
+    /*  might hold one of these in a register. */
+
+    /*  Arguments */
+    /*  ========= */
+
+    /*  A       (input) DOUBLE PRECISION */
+    /*  B       (input) DOUBLE PRECISION */
+    /*          The values A and B. */
+
+    /* ===================================================================== */
+
+    /*     .. Executable Statements .. */
+
+    double ret_val = *a + *b;
+
+    return ret_val;
+
+    /*     End of DLAMC3 */
+
+} /* dlamc3_ */
+
+
+/* simpler version of dlamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
+   taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
+
+#ifndef DBL_DIGITS
+#define DBL_DIGITS 53
+#endif
+
+static const unsigned char lapack_dlamch_tab0[] =
+{
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9, 0, 10, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9,
+    0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+};
+
+const double lapack_dlamch_tab1[] =
+{
+    0, FLT_RADIX, DBL_EPSILON, DBL_MAX_EXP, DBL_MIN_EXP, DBL_DIGITS, DBL_MAX,
+    DBL_EPSILON*FLT_RADIX, 1, DBL_MIN*(1 + DBL_EPSILON), DBL_MIN
+};
+
+double dlamch_(char* cmach)
+{
+    return lapack_dlamch_tab1[lapack_dlamch_tab0[(unsigned char)cmach[0]]];
+}
@@ -0,0 +1,96 @@
+#include "f2c.h"
+
+static const int CLAPACK_NOT_IMPLEMENTED = -1024;
+
+int sgesdd_(char *jobz, int *m, int *n, float *a, int *lda,
+	float *s, float *u, int *ldu, float *vt, int *ldvt, float *work,
+	int *lwork, int *iwork, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int dgels_(char *trans, int *m, int *n, int *nrhs, double *a,
+	int *lda, double *b, int *ldb, double *work, int *lwork, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv,
+	double *b, int *ldb, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv,
+	int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int dposv_(char *uplo, int *n, int *nrhs, double *a, int *
+	lda, double *b, int *ldb, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int dpotrf_(char *uplo, int *n, double *a, int *lda, int *
+	info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int sgels_(char *trans, int *m, int *n, int *nrhs, float *a,
+	int *lda, float *b, int *ldb, float *work, int *lwork, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int sgeev_(char *jobvl, char *jobvr, int *n, float *a, int *
+	lda, float *wr, float *wi, float *vl, int *ldvl, float *vr, int *
+	ldvr, float *work, int *lwork, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int sgeqrf_(int *m, int *n, float *a, int *lda, float *tau,
+	float *work, int *lwork, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int sgesv_(int *n, int *nrhs, float *a, int *lda, int *ipiv,
+	float *b, int *ldb, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}        
+
+int sgetrf_(int *m, int *n, float *a, int *lda, int *ipiv,
+	int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
+
+int sposv_(char *uplo, int *n, int *nrhs, float *a, int *
+	lda, float *b, int *ldb, int *info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}        
+
+int spotrf_(char *uplo, int *n, float *a, int *lda, int *
+	info)
+{
+    *info = CLAPACK_NOT_IMPLEMENTED;
+    return 0;
+}
@@ -0,0 +1,25 @@
+#include "f2c.h"
+
+static const unsigned char lapack_toupper_tab[] =
+{
+    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
+    24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
+    46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
+    68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
+    90, 91, 92, 93, 94, 95, 96, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
+    80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 123, 124, 125, 126, 127, 128, 129, 130, 131,
+    132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
+    150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167,
+    168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185,
+    186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203,
+    204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221,
+    222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
+    240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255
+};
+
+#define lapack_toupper(c) ((char)lapack_toupper_tab[(unsigned char)(c)])
+
+int lsame_(char *ca, char *cb)
+{
+    return lapack_toupper(ca[0]) == lapack_toupper(cb[0]);
+}
@@ -0,0 +1,27 @@
+#include "f2c.h"
+
+double pow_di(double *ap, int *bp)
+{
+    double p = 1;
+    double x = *ap;
+    int n = *bp;
+
+    if(n != 0)
+    {
+        if(n < 0)
+        {
+            n = -n;
+            x = 1/x;
+        }
+        unsigned u = (unsigned)n;
+        for(;;)
+        {
+            if((u & 1) != 0)
+                p *= x;
+            if((u >>= 1) == 0)
+                break;
+            x *= x;
+        }
+    }
+    return p;
+}
@@ -0,0 +1,25 @@
+#include "f2c.h"
+
+int pow_ii(int *ap, int *bp)
+{
+    int p;
+    int x = *ap;
+    int n = *bp;
+
+    if (n <= 0) {
+        if (n == 0 || x == 1)
+            return 1;
+        return x != -1 ? 0 : (n & 1) ? -1 : 1;
+    }
+    unsigned u = (unsigned)n;
+    for(p = 1; ; )
+    {
+        if(u & 01)
+            p *= x;
+        if(u >>= 1)
+            x *= x;
+        else
+            break;
+    }
+    return p;
+}
@@ -0,0 +1,22 @@
+/* Unless compiled with -DNO_OVERWRITE, this variant of s_cat allows the
+ * target of a concatenation to appear on its right-hand side (contrary
+ * to the Fortran 77 Standard, but in accordance with Fortran 90).
+ */
+
+#include "f2c.h"
+
+int s_cat(char *lp, char **rpp, int* rnp, int *np)
+{
+    int i, L = 0;
+    int n = *np;
+
+    for(i = 0; i < n; i++) {
+        int ni = rnp[i];
+        if(ni > 0) {
+            memcpy(lp + L, rpp[i], ni);
+            L += ni;
+        }
+    }
+    lp[L] = '\0';
+    return 0;
+}
@@ -0,0 +1,40 @@
+#include "f2c.h"
+
+/* compare two strings */
+int s_cmp(char *a0, char *b0)
+{
+    int la = (int)strlen(a0);
+    int lb = (int)strlen(b0);
+    unsigned char *a, *aend, *b, *bend;
+    a = (unsigned char *)a0;
+    b = (unsigned char *)b0;
+    aend = a + la;
+    bend = b + lb;
+
+    if(la <= lb)
+    	{
+    	while(a < aend)
+    		if(*a != *b)
+    			return( *a - *b );
+    		else
+    			{ ++a; ++b; }
+
+    	while(b < bend)
+    		if(*b != ' ')
+    			return( ' ' - *b );
+    		else	++b;
+    	}
+    else
+    	{
+    	while(b < bend)
+    		if(*a == *b)
+    			{ ++a; ++b; }
+    		else
+    			return( *a - *b );
+    	while(a < aend)
+    		if(*a != ' ')
+    			return(*a - ' ');
+    		else	++a;
+    	}
+    return(0);
+}
@@ -0,0 +1,71 @@
+#include "f2c.h"
+#include <float.h>
+#include <stdio.h>
+
+/* *********************************************************************** */
+
+double slamc3_(float *a, float *b)
+{
+    /*  -- LAPACK auxiliary routine (version 3.1) -- */
+    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+    /*     November 2006 */
+
+    /*     .. Scalar Arguments .. */
+    /*     .. */
+
+    /*  Purpose */
+    /*  ======= */
+
+    /*  SLAMC3  is intended to force  A  and  B  to be stored prior to doing */
+    /*  the addition of  A  and  B ,  for use in situations where optimizers */
+    /*  might hold one of these in a register. */
+
+    /*  Arguments */
+    /*  ========= */
+
+    /*  A       (input) REAL */
+    /*  B       (input) REAL */
+    /*          The values A and B. */
+
+    /* ===================================================================== */
+
+    /*     .. Executable Statements .. */
+
+    float ret_val = *a + *b;
+
+    return ret_val;
+
+    /*     End of SLAMC3 */
+
+} /* slamc3_ */
+
+/* simpler version of slamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
+ taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
+
+#ifndef FLT_DIGITS
+#define FLT_DIGITS 24
+#endif
+
+static const unsigned char lapack_slamch_tab0[] =
+{
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9, 0, 10, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 4, 5, 6, 7, 0, 8, 9,
+    0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+};
+
+const double lapack_slamch_tab1[] =
+{
+    0, FLT_RADIX, FLT_EPSILON, FLT_MAX_EXP, FLT_MIN_EXP, FLT_DIGITS, FLT_MAX,
+    FLT_EPSILON*FLT_RADIX, 1, FLT_MIN*(1 + FLT_EPSILON), FLT_MIN
+};
+
+double slamch_(char* cmach)
+{
+    return lapack_slamch_tab1[lapack_slamch_tab0[(unsigned char)cmach[0]]];
+}
@@ -0,0 +1,19 @@
+/* xerbla.f -- translated by f2c (version 20061008).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "f2c.h"
+
+/* Subroutine */ int xerbla_(char *srname, int *info)
+{
+    printf("** On entry to %s, parameter number %2i had an illegal value\n", srname, *info);
+    return 0;
+} /* xerbla_ */
@@ -0,0 +1,752 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b CGEMM
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE CGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+//
+//      .. Scalar Arguments ..
+//      COMPLEX ALPHA,BETA
+//      INTEGER K,LDA,LDB,LDC,M,N
+//      CHARACTER TRANSA,TRANSB
+//      ..
+//      .. Array Arguments ..
+//      COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> CGEMM  performs one of the matrix-matrix operations
+//>
+//>    C := alpha*op( A )*op( B ) + beta*C,
+//>
+//> where  op( X ) is one of
+//>
+//>    op( X ) = X   or   op( X ) = X**T   or   op( X ) = X**H,
+//>
+//> alpha and beta are scalars, and A, B and C are matrices, with op( A )
+//> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] TRANSA
+//> \verbatim
+//>          TRANSA is CHARACTER*1
+//>           On entry, TRANSA specifies the form of op( A ) to be used in
+//>           the matrix multiplication as follows:
+//>
+//>              TRANSA = 'N' or 'n',  op( A ) = A.
+//>
+//>              TRANSA = 'T' or 't',  op( A ) = A**T.
+//>
+//>              TRANSA = 'C' or 'c',  op( A ) = A**H.
+//> \endverbatim
+//>
+//> \param[in] TRANSB
+//> \verbatim
+//>          TRANSB is CHARACTER*1
+//>           On entry, TRANSB specifies the form of op( B ) to be used in
+//>           the matrix multiplication as follows:
+//>
+//>              TRANSB = 'N' or 'n',  op( B ) = B.
+//>
+//>              TRANSB = 'T' or 't',  op( B ) = B**T.
+//>
+//>              TRANSB = 'C' or 'c',  op( B ) = B**H.
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>           On entry,  M  specifies  the number  of rows  of the  matrix
+//>           op( A )  and of the  matrix  C.  M  must  be at least  zero.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>           On entry,  N  specifies the number  of columns of the matrix
+//>           op( B ) and the number of columns of the matrix C. N must be
+//>           at least zero.
+//> \endverbatim
+//>
+//> \param[in] K
+//> \verbatim
+//>          K is INTEGER
+//>           On entry,  K  specifies  the number of columns of the matrix
+//>           op( A ) and the number of rows of the matrix op( B ). K must
+//>           be at least  zero.
+//> \endverbatim
+//>
+//> \param[in] ALPHA
+//> \verbatim
+//>          ALPHA is COMPLEX
+//>           On entry, ALPHA specifies the scalar alpha.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is COMPLEX array, dimension ( LDA, ka ), where ka is
+//>           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
+//>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
+//>           part of the array  A  must contain the matrix  A,  otherwise
+//>           the leading  k by m  part of the array  A  must contain  the
+//>           matrix A.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>           On entry, LDA specifies the first dimension of A as declared
+//>           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
+//>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
+//>           least  max( 1, k ).
+//> \endverbatim
+//>
+//> \param[in] B
+//> \verbatim
+//>          B is COMPLEX array, dimension ( LDB, kb ), where kb is
+//>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
+//>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
+//>           part of the array  B  must contain the matrix  B,  otherwise
+//>           the leading  n by k  part of the array  B  must contain  the
+//>           matrix B.
+//> \endverbatim
+//>
+//> \param[in] LDB
+//> \verbatim
+//>          LDB is INTEGER
+//>           On entry, LDB specifies the first dimension of B as declared
+//>           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
+//>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
+//>           least  max( 1, n ).
+//> \endverbatim
+//>
+//> \param[in] BETA
+//> \verbatim
+//>          BETA is COMPLEX
+//>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
+//>           supplied as zero then C need not be set on input.
+//> \endverbatim
+//>
+//> \param[in,out] C
+//> \verbatim
+//>          C is COMPLEX array, dimension ( LDC, N )
+//>           Before entry, the leading  m by n  part of the array  C must
+//>           contain the matrix  C,  except when  beta  is zero, in which
+//>           case C need not be set on entry.
+//>           On exit, the array  C  is overwritten by the  m by n  matrix
+//>           ( alpha*op( A )*op( B ) + beta*C ).
+//> \endverbatim
+//>
+//> \param[in] LDC
+//> \verbatim
+//>          LDC is INTEGER
+//>           On entry, LDC specifies the first dimension of C as declared
+//>           in  the  calling  (sub)  program.   LDC  must  be  at  least
+//>           max( 1, m ).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup complex_blas_level3
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>  Level 3 Blas routine.
+//>
+//>  -- Written on 8-February-1989.
+//>     Jack Dongarra, Argonne National Laboratory.
+//>     Iain Duff, AERE Harwell.
+//>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
+//>     Sven Hammarling, Numerical Algorithms Group Ltd.
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int cgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, complex *alpha, complex *a, int *lda, complex *b, int *ldb,
+	complex *beta, complex *c__, int *ldc)
+{
+    // Table of constant values
+    complex c_b1 = {1.f,0.f};
+    complex c_b2 = {0.f,0.f};
+
+    // System generated locals
+    int a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
+	    i__3, i__4, i__5, i__6;
+    complex q__1, q__2, q__3, q__4;
+
+    // Local variables
+    int i__, j, l, info;
+    int nota, notb;
+    complex temp;
+    int conja, conjb;
+    int ncola;
+    extern int lsame_(char *, char *);
+    int nrowa, nrowb;
+    extern /* Subroutine */ int xerbla_(char *, int *);
+
+    //
+    // -- Reference BLAS level3 routine (version 3.7.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. External Functions ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Parameters ..
+    //    ..
+    //
+    //    Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
+    //    conjugated or transposed, set  CONJA and CONJB  as true if  A  and
+    //    B  respectively are to be  transposed but  not conjugated  and set
+    //    NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A
+    //    and the number of rows of  B  respectively.
+    //
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    b_dim1 = *ldb;
+    b_offset = 1 + b_dim1;
+    b -= b_offset;
+    c_dim1 = *ldc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    // Function Body
+    nota = lsame_(transa, "N");
+    notb = lsame_(transb, "N");
+    conja = lsame_(transa, "C");
+    conjb = lsame_(transb, "C");
+    if (nota) {
+	nrowa = *m;
+	ncola = *k;
+    } else {
+	nrowa = *k;
+	ncola = *m;
+    }
+    if (notb) {
+	nrowb = *k;
+    } else {
+	nrowb = *n;
+    }
+    //
+    //    Test the input parameters.
+    //
+    info = 0;
+    if (! nota && ! conja && ! lsame_(transa, "T")) {
+	info = 1;
+    } else if (! notb && ! conjb && ! lsame_(transb, "T")) {
+	info = 2;
+    } else if (*m < 0) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < max(1,nrowa)) {
+	info = 8;
+    } else if (*ldb < max(1,nrowb)) {
+	info = 10;
+    } else if (*ldc < max(1,*m)) {
+	info = 13;
+    }
+    if (info != 0) {
+	xerbla_("CGEMM ", &info);
+	return 0;
+    }
+    //
+    //    Quick return if possible.
+    //
+    if (*m == 0 || *n == 0 || (alpha->r == 0.f && alpha->i == 0.f || *k == 0)
+	    && (beta->r == 1.f && beta->i == 0.f)) {
+	return 0;
+    }
+    //
+    //    And when  alpha.eq.zero.
+    //
+    if (alpha->r == 0.f && alpha->i == 0.f) {
+	if (beta->r == 0.f && beta->i == 0.f) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    i__3 = i__ + j * c_dim1;
+		    c__[i__3].r = 0.f, c__[i__3].i = 0.f;
+// L10:
+		}
+// L20:
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    i__3 = i__ + j * c_dim1;
+		    i__4 = i__ + j * c_dim1;
+		    q__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4].i,
+			    q__1.i = beta->r * c__[i__4].i + beta->i * c__[
+			    i__4].r;
+		    c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+// L30:
+		}
+// L40:
+	    }
+	}
+	return 0;
+    }
+    //
+    //    Start the operations.
+    //
+    if (notb) {
+	if (nota) {
+	    //
+	    //          Form  C := alpha*A*B + beta*C.
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (beta->r == 0.f && beta->i == 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			c__[i__3].r = 0.f, c__[i__3].i = 0.f;
+// L50:
+		    }
+		} else if (beta->r != 1.f || beta->i != 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			i__4 = i__ + j * c_dim1;
+			q__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__1.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+// L60:
+		    }
+		}
+		i__2 = *k;
+		for (l = 1; l <= i__2; ++l) {
+		    i__3 = l + j * b_dim1;
+		    q__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
+			    q__1.i = alpha->r * b[i__3].i + alpha->i * b[i__3]
+			    .r;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    i__3 = *m;
+		    for (i__ = 1; i__ <= i__3; ++i__) {
+			i__4 = i__ + j * c_dim1;
+			i__5 = i__ + j * c_dim1;
+			i__6 = i__ + l * a_dim1;
+			q__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
+				q__2.i = temp.r * a[i__6].i + temp.i * a[i__6]
+				.r;
+			q__1.r = c__[i__5].r + q__2.r, q__1.i = c__[i__5].i +
+				q__2.i;
+			c__[i__4].r = q__1.r, c__[i__4].i = q__1.i;
+// L70:
+		    }
+// L80:
+		}
+// L90:
+	    }
+	} else if (conja) {
+	    //
+	    //          Form  C := alpha*A**H*B + beta*C.
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			r_cnjg(&q__3, &a[l + i__ * a_dim1]);
+			i__4 = l + j * b_dim1;
+			q__2.r = q__3.r * b[i__4].r - q__3.i * b[i__4].i,
+				q__2.i = q__3.r * b[i__4].i + q__3.i * b[i__4]
+				.r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L100:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L110:
+		}
+// L120:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A**T*B + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			i__4 = l + i__ * a_dim1;
+			i__5 = l + j * b_dim1;
+			q__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
+				.i, q__2.i = a[i__4].r * b[i__5].i + a[i__4]
+				.i * b[i__5].r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L130:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L140:
+		}
+// L150:
+	    }
+	}
+    } else if (nota) {
+	if (conjb) {
+	    //
+	    //          Form  C := alpha*A*B**H + beta*C.
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (beta->r == 0.f && beta->i == 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			c__[i__3].r = 0.f, c__[i__3].i = 0.f;
+// L160:
+		    }
+		} else if (beta->r != 1.f || beta->i != 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			i__4 = i__ + j * c_dim1;
+			q__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__1.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+// L170:
+		    }
+		}
+		i__2 = *k;
+		for (l = 1; l <= i__2; ++l) {
+		    r_cnjg(&q__2, &b[j + l * b_dim1]);
+		    q__1.r = alpha->r * q__2.r - alpha->i * q__2.i, q__1.i =
+			    alpha->r * q__2.i + alpha->i * q__2.r;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    i__3 = *m;
+		    for (i__ = 1; i__ <= i__3; ++i__) {
+			i__4 = i__ + j * c_dim1;
+			i__5 = i__ + j * c_dim1;
+			i__6 = i__ + l * a_dim1;
+			q__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
+				q__2.i = temp.r * a[i__6].i + temp.i * a[i__6]
+				.r;
+			q__1.r = c__[i__5].r + q__2.r, q__1.i = c__[i__5].i +
+				q__2.i;
+			c__[i__4].r = q__1.r, c__[i__4].i = q__1.i;
+// L180:
+		    }
+// L190:
+		}
+// L200:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A*B**T + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (beta->r == 0.f && beta->i == 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			c__[i__3].r = 0.f, c__[i__3].i = 0.f;
+// L210:
+		    }
+		} else if (beta->r != 1.f || beta->i != 0.f) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			i__3 = i__ + j * c_dim1;
+			i__4 = i__ + j * c_dim1;
+			q__1.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__1.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+// L220:
+		    }
+		}
+		i__2 = *k;
+		for (l = 1; l <= i__2; ++l) {
+		    i__3 = j + l * b_dim1;
+		    q__1.r = alpha->r * b[i__3].r - alpha->i * b[i__3].i,
+			    q__1.i = alpha->r * b[i__3].i + alpha->i * b[i__3]
+			    .r;
+		    temp.r = q__1.r, temp.i = q__1.i;
+		    i__3 = *m;
+		    for (i__ = 1; i__ <= i__3; ++i__) {
+			i__4 = i__ + j * c_dim1;
+			i__5 = i__ + j * c_dim1;
+			i__6 = i__ + l * a_dim1;
+			q__2.r = temp.r * a[i__6].r - temp.i * a[i__6].i,
+				q__2.i = temp.r * a[i__6].i + temp.i * a[i__6]
+				.r;
+			q__1.r = c__[i__5].r + q__2.r, q__1.i = c__[i__5].i +
+				q__2.i;
+			c__[i__4].r = q__1.r, c__[i__4].i = q__1.i;
+// L230:
+		    }
+// L240:
+		}
+// L250:
+	    }
+	}
+    } else if (conja) {
+	if (conjb) {
+	    //
+	    //          Form  C := alpha*A**H*B**H + beta*C.
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			r_cnjg(&q__3, &a[l + i__ * a_dim1]);
+			r_cnjg(&q__4, &b[j + l * b_dim1]);
+			q__2.r = q__3.r * q__4.r - q__3.i * q__4.i, q__2.i =
+				q__3.r * q__4.i + q__3.i * q__4.r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L260:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L270:
+		}
+// L280:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A**H*B**T + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			r_cnjg(&q__3, &a[l + i__ * a_dim1]);
+			i__4 = j + l * b_dim1;
+			q__2.r = q__3.r * b[i__4].r - q__3.i * b[i__4].i,
+				q__2.i = q__3.r * b[i__4].i + q__3.i * b[i__4]
+				.r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L290:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L300:
+		}
+// L310:
+	    }
+	}
+    } else {
+	if (conjb) {
+	    //
+	    //          Form  C := alpha*A**T*B**H + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			i__4 = l + i__ * a_dim1;
+			r_cnjg(&q__3, &b[j + l * b_dim1]);
+			q__2.r = a[i__4].r * q__3.r - a[i__4].i * q__3.i,
+				q__2.i = a[i__4].r * q__3.i + a[i__4].i *
+				q__3.r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L320:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L330:
+		}
+// L340:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A**T*B**T + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp.r = 0.f, temp.i = 0.f;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			i__4 = l + i__ * a_dim1;
+			i__5 = j + l * b_dim1;
+			q__2.r = a[i__4].r * b[i__5].r - a[i__4].i * b[i__5]
+				.i, q__2.i = a[i__4].r * b[i__5].i + a[i__4]
+				.i * b[i__5].r;
+			q__1.r = temp.r + q__2.r, q__1.i = temp.i + q__2.i;
+			temp.r = q__1.r, temp.i = q__1.i;
+// L350:
+		    }
+		    if (beta->r == 0.f && beta->i == 0.f) {
+			i__3 = i__ + j * c_dim1;
+			q__1.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__1.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    } else {
+			i__3 = i__ + j * c_dim1;
+			q__2.r = alpha->r * temp.r - alpha->i * temp.i,
+				q__2.i = alpha->r * temp.i + alpha->i *
+				temp.r;
+			i__4 = i__ + j * c_dim1;
+			q__3.r = beta->r * c__[i__4].r - beta->i * c__[i__4]
+				.i, q__3.i = beta->r * c__[i__4].i + beta->i *
+				 c__[i__4].r;
+			q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
+			c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
+		    }
+// L360:
+		}
+// L370:
+	    }
+	}
+    }
+    return 0;
+    //
+    //    End of CGEMM .
+    //
+} // cgemm_
+
@@ -0,0 +1,171 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DCOPY
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
+//
+//      .. Scalar Arguments ..
+//      INTEGER INCX,INCY,N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION DX(*),DY(*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//>    DCOPY copies a vector, x, to a vector, y.
+//>    uses unrolled loops for increments equal to 1.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>         number of elements in input vector(s)
+//> \endverbatim
+//>
+//> \param[in] DX
+//> \verbatim
+//>          DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+//> \endverbatim
+//>
+//> \param[in] INCX
+//> \verbatim
+//>          INCX is INTEGER
+//>         storage spacing between elements of DX
+//> \endverbatim
+//>
+//> \param[out] DY
+//> \verbatim
+//>          DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+//> \endverbatim
+//>
+//> \param[in] INCY
+//> \verbatim
+//>          INCY is INTEGER
+//>         storage spacing between elements of DY
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date November 2017
+//
+//> \ingroup double_blas_level1
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>     jack dongarra, linpack, 3/11/78.
+//>     modified 12/3/93, array(1) declarations changed to array(*)
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int dcopy_(int *n, double *dx, int *incx, double *dy, int *
+	incy)
+{
+    // System generated locals
+    int i__1;
+
+    // Local variables
+    int i__, m, ix, iy, mp1;
+
+    //
+    // -- Reference BLAS level1 routine (version 3.8.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    November 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    // Parameter adjustments
+    --dy;
+    --dx;
+
+    // Function Body
+    if (*n <= 0) {
+	return 0;
+    }
+    if (*incx == 1 && *incy == 1) {
+	//
+	//       code for both increments equal to 1
+	//
+	//
+	//       clean-up loop
+	//
+	m = *n % 7;
+	if (m != 0) {
+	    i__1 = m;
+	    for (i__ = 1; i__ <= i__1; ++i__) {
+		dy[i__] = dx[i__];
+	    }
+	    if (*n < 7) {
+		return 0;
+	    }
+	}
+	mp1 = m + 1;
+	i__1 = *n;
+	for (i__ = mp1; i__ <= i__1; i__ += 7) {
+	    dy[i__] = dx[i__];
+	    dy[i__ + 1] = dx[i__ + 1];
+	    dy[i__ + 2] = dx[i__ + 2];
+	    dy[i__ + 3] = dx[i__ + 3];
+	    dy[i__ + 4] = dx[i__ + 4];
+	    dy[i__ + 5] = dx[i__ + 5];
+	    dy[i__ + 6] = dx[i__ + 6];
+	}
+    } else {
+	//
+	//       code for unequal increments or equal increments
+	//         not equal to 1
+	//
+	ix = 1;
+	iy = 1;
+	if (*incx < 0) {
+	    ix = (-(*n) + 1) * *incx + 1;
+	}
+	if (*incy < 0) {
+	    iy = (-(*n) + 1) * *incy + 1;
+	}
+	i__1 = *n;
+	for (i__ = 1; i__ <= i__1; ++i__) {
+	    dy[iy] = dx[ix];
+	    ix += *incx;
+	    iy += *incy;
+	}
+    }
+    return 0;
+} // dcopy_
+
@@ -0,0 +1,172 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DDOT
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
+//
+//      .. Scalar Arguments ..
+//      INTEGER INCX,INCY,N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION DX(*),DY(*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//>    DDOT forms the dot product of two vectors.
+//>    uses unrolled loops for increments equal to one.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>         number of elements in input vector(s)
+//> \endverbatim
+//>
+//> \param[in] DX
+//> \verbatim
+//>          DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+//> \endverbatim
+//>
+//> \param[in] INCX
+//> \verbatim
+//>          INCX is INTEGER
+//>         storage spacing between elements of DX
+//> \endverbatim
+//>
+//> \param[in] DY
+//> \verbatim
+//>          DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+//> \endverbatim
+//>
+//> \param[in] INCY
+//> \verbatim
+//>          INCY is INTEGER
+//>         storage spacing between elements of DY
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date November 2017
+//
+//> \ingroup double_blas_level1
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>     jack dongarra, linpack, 3/11/78.
+//>     modified 12/3/93, array(1) declarations changed to array(*)
+//> \endverbatim
+//>
+// =====================================================================
+double ddot_(int *n, double *dx, int *incx, double *dy, int *incy)
+{
+    // System generated locals
+    int i__1;
+    double ret_val;
+
+    // Local variables
+    int i__, m, ix, iy, mp1;
+    double dtemp;
+
+    //
+    // -- Reference BLAS level1 routine (version 3.8.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    November 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    // Parameter adjustments
+    --dy;
+    --dx;
+
+    // Function Body
+    ret_val = 0.;
+    dtemp = 0.;
+    if (*n <= 0) {
+	return ret_val;
+    }
+    if (*incx == 1 && *incy == 1) {
+	//
+	//       code for both increments equal to 1
+	//
+	//
+	//       clean-up loop
+	//
+	m = *n % 5;
+	if (m != 0) {
+	    i__1 = m;
+	    for (i__ = 1; i__ <= i__1; ++i__) {
+		dtemp += dx[i__] * dy[i__];
+	    }
+	    if (*n < 5) {
+		ret_val = dtemp;
+		return ret_val;
+	    }
+	}
+	mp1 = m + 1;
+	i__1 = *n;
+	for (i__ = mp1; i__ <= i__1; i__ += 5) {
+	    dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] +
+		    dx[i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] +
+		    dx[i__ + 4] * dy[i__ + 4];
+	}
+    } else {
+	//
+	//       code for unequal increments or equal increments
+	//         not equal to 1
+	//
+	ix = 1;
+	iy = 1;
+	if (*incx < 0) {
+	    ix = (-(*n) + 1) * *incx + 1;
+	}
+	if (*incy < 0) {
+	    iy = (-(*n) + 1) * *incy + 1;
+	}
+	i__1 = *n;
+	for (i__ = 1; i__ <= i__1; ++i__) {
+	    dtemp += dx[ix] * dy[iy];
+	    ix += *incx;
+	    iy += *incy;
+	}
+    }
+    ret_val = dtemp;
+    return ret_val;
+} // ddot_
+
@@ -0,0 +1,444 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DGEMM
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION ALPHA,BETA
+//      INTEGER K,LDA,LDB,LDC,M,N
+//      CHARACTER TRANSA,TRANSB
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DGEMM  performs one of the matrix-matrix operations
+//>
+//>    C := alpha*op( A )*op( B ) + beta*C,
+//>
+//> where  op( X ) is one of
+//>
+//>    op( X ) = X   or   op( X ) = X**T,
+//>
+//> alpha and beta are scalars, and A, B and C are matrices, with op( A )
+//> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] TRANSA
+//> \verbatim
+//>          TRANSA is CHARACTER*1
+//>           On entry, TRANSA specifies the form of op( A ) to be used in
+//>           the matrix multiplication as follows:
+//>
+//>              TRANSA = 'N' or 'n',  op( A ) = A.
+//>
+//>              TRANSA = 'T' or 't',  op( A ) = A**T.
+//>
+//>              TRANSA = 'C' or 'c',  op( A ) = A**T.
+//> \endverbatim
+//>
+//> \param[in] TRANSB
+//> \verbatim
+//>          TRANSB is CHARACTER*1
+//>           On entry, TRANSB specifies the form of op( B ) to be used in
+//>           the matrix multiplication as follows:
+//>
+//>              TRANSB = 'N' or 'n',  op( B ) = B.
+//>
+//>              TRANSB = 'T' or 't',  op( B ) = B**T.
+//>
+//>              TRANSB = 'C' or 'c',  op( B ) = B**T.
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>           On entry,  M  specifies  the number  of rows  of the  matrix
+//>           op( A )  and of the  matrix  C.  M  must  be at least  zero.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>           On entry,  N  specifies the number  of columns of the matrix
+//>           op( B ) and the number of columns of the matrix C. N must be
+//>           at least zero.
+//> \endverbatim
+//>
+//> \param[in] K
+//> \verbatim
+//>          K is INTEGER
+//>           On entry,  K  specifies  the number of columns of the matrix
+//>           op( A ) and the number of rows of the matrix op( B ). K must
+//>           be at least  zero.
+//> \endverbatim
+//>
+//> \param[in] ALPHA
+//> \verbatim
+//>          ALPHA is DOUBLE PRECISION.
+//>           On entry, ALPHA specifies the scalar alpha.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
+//>           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
+//>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
+//>           part of the array  A  must contain the matrix  A,  otherwise
+//>           the leading  k by m  part of the array  A  must contain  the
+//>           matrix A.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>           On entry, LDA specifies the first dimension of A as declared
+//>           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
+//>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
+//>           least  max( 1, k ).
+//> \endverbatim
+//>
+//> \param[in] B
+//> \verbatim
+//>          B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
+//>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
+//>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
+//>           part of the array  B  must contain the matrix  B,  otherwise
+//>           the leading  n by k  part of the array  B  must contain  the
+//>           matrix B.
+//> \endverbatim
+//>
+//> \param[in] LDB
+//> \verbatim
+//>          LDB is INTEGER
+//>           On entry, LDB specifies the first dimension of B as declared
+//>           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
+//>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
+//>           least  max( 1, n ).
+//> \endverbatim
+//>
+//> \param[in] BETA
+//> \verbatim
+//>          BETA is DOUBLE PRECISION.
+//>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
+//>           supplied as zero then C need not be set on input.
+//> \endverbatim
+//>
+//> \param[in,out] C
+//> \verbatim
+//>          C is DOUBLE PRECISION array, dimension ( LDC, N )
+//>           Before entry, the leading  m by n  part of the array  C must
+//>           contain the matrix  C,  except when  beta  is zero, in which
+//>           case C need not be set on entry.
+//>           On exit, the array  C  is overwritten by the  m by n  matrix
+//>           ( alpha*op( A )*op( B ) + beta*C ).
+//> \endverbatim
+//>
+//> \param[in] LDC
+//> \verbatim
+//>          LDC is INTEGER
+//>           On entry, LDC specifies the first dimension of C as declared
+//>           in  the  calling  (sub)  program.   LDC  must  be  at  least
+//>           max( 1, m ).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup double_blas_level3
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>  Level 3 Blas routine.
+//>
+//>  -- Written on 8-February-1989.
+//>     Jack Dongarra, Argonne National Laboratory.
+//>     Iain Duff, AERE Harwell.
+//>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
+//>     Sven Hammarling, Numerical Algorithms Group Ltd.
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int dgemm_(char *transa, char *transb, int *m, int *n, int *
+	k, double *alpha, double *a, int *lda, double *b, int *ldb, double *
+	beta, double *c__, int *ldc)
+{
+    // System generated locals
+    int a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
+	    i__3;
+
+    // Local variables
+    int i__, j, l, info;
+    int nota, notb;
+    double temp;
+    int ncola;
+    extern int lsame_(char *, char *);
+    int nrowa, nrowb;
+    extern /* Subroutine */ int xerbla_(char *, int *);
+
+    //
+    // -- Reference BLAS level3 routine (version 3.7.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. External Functions ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Parameters ..
+    //    ..
+    //
+    //    Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
+    //    transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
+    //    and  columns of  A  and the  number of  rows  of  B  respectively.
+    //
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    b_dim1 = *ldb;
+    b_offset = 1 + b_dim1;
+    b -= b_offset;
+    c_dim1 = *ldc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    // Function Body
+    nota = lsame_(transa, "N");
+    notb = lsame_(transb, "N");
+    if (nota) {
+	nrowa = *m;
+	ncola = *k;
+    } else {
+	nrowa = *k;
+	ncola = *m;
+    }
+    if (notb) {
+	nrowb = *k;
+    } else {
+	nrowb = *n;
+    }
+    //
+    //    Test the input parameters.
+    //
+    info = 0;
+    if (! nota && ! lsame_(transa, "C") && ! lsame_(transa, "T")) {
+	info = 1;
+    } else if (! notb && ! lsame_(transb, "C") && ! lsame_(transb, "T")) {
+	info = 2;
+    } else if (*m < 0) {
+	info = 3;
+    } else if (*n < 0) {
+	info = 4;
+    } else if (*k < 0) {
+	info = 5;
+    } else if (*lda < max(1,nrowa)) {
+	info = 8;
+    } else if (*ldb < max(1,nrowb)) {
+	info = 10;
+    } else if (*ldc < max(1,*m)) {
+	info = 13;
+    }
+    if (info != 0) {
+	xerbla_("DGEMM ", &info);
+	return 0;
+    }
+    //
+    //    Quick return if possible.
+    //
+    if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
+	return 0;
+    }
+    //
+    //    And if  alpha.eq.zero.
+    //
+    if (*alpha == 0.) {
+	if (*beta == 0.) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    c__[i__ + j * c_dim1] = 0.;
+// L10:
+		}
+// L20:
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
+// L30:
+		}
+// L40:
+	    }
+	}
+	return 0;
+    }
+    //
+    //    Start the operations.
+    //
+    if (notb) {
+	if (nota) {
+	    //
+	    //          Form  C := alpha*A*B + beta*C.
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (*beta == 0.) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] = 0.;
+// L50:
+		    }
+		} else if (*beta != 1.) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
+// L60:
+		    }
+		}
+		i__2 = *k;
+		for (l = 1; l <= i__2; ++l) {
+		    temp = *alpha * b[l + j * b_dim1];
+		    i__3 = *m;
+		    for (i__ = 1; i__ <= i__3; ++i__) {
+			c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1];
+// L70:
+		    }
+// L80:
+		}
+// L90:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A**T*B + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp = 0.;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
+// L100:
+		    }
+		    if (*beta == 0.) {
+			c__[i__ + j * c_dim1] = *alpha * temp;
+		    } else {
+			c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
+				i__ + j * c_dim1];
+		    }
+// L110:
+		}
+// L120:
+	    }
+	}
+    } else {
+	if (nota) {
+	    //
+	    //          Form  C := alpha*A*B**T + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		if (*beta == 0.) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] = 0.;
+// L130:
+		    }
+		} else if (*beta != 1.) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
+// L140:
+		    }
+		}
+		i__2 = *k;
+		for (l = 1; l <= i__2; ++l) {
+		    temp = *alpha * b[j + l * b_dim1];
+		    i__3 = *m;
+		    for (i__ = 1; i__ <= i__3; ++i__) {
+			c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1];
+// L150:
+		    }
+// L160:
+		}
+// L170:
+	    }
+	} else {
+	    //
+	    //          Form  C := alpha*A**T*B**T + beta*C
+	    //
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp = 0.;
+		    i__3 = *k;
+		    for (l = 1; l <= i__3; ++l) {
+			temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
+// L180:
+		    }
+		    if (*beta == 0.) {
+			c__[i__ + j * c_dim1] = *alpha * temp;
+		    } else {
+			c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
+				i__ + j * c_dim1];
+		    }
+// L190:
+		}
+// L200:
+	    }
+	}
+    }
+    return 0;
+    //
+    //    End of DGEMM .
+    //
+} // dgemm_
+
@@ -0,0 +1,370 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DGEMV
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION ALPHA,BETA
+//      INTEGER INCX,INCY,LDA,M,N
+//      CHARACTER TRANS
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DGEMV  performs one of the matrix-vector operations
+//>
+//>    y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,
+//>
+//> where alpha and beta are scalars, x and y are vectors and A is an
+//> m by n matrix.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] TRANS
+//> \verbatim
+//>          TRANS is CHARACTER*1
+//>           On entry, TRANS specifies the operation to be performed as
+//>           follows:
+//>
+//>              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
+//>
+//>              TRANS = 'T' or 't'   y := alpha*A**T*x + beta*y.
+//>
+//>              TRANS = 'C' or 'c'   y := alpha*A**T*x + beta*y.
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>           On entry, M specifies the number of rows of the matrix A.
+//>           M must be at least zero.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>           On entry, N specifies the number of columns of the matrix A.
+//>           N must be at least zero.
+//> \endverbatim
+//>
+//> \param[in] ALPHA
+//> \verbatim
+//>          ALPHA is DOUBLE PRECISION.
+//>           On entry, ALPHA specifies the scalar alpha.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension ( LDA, N )
+//>           Before entry, the leading m by n part of the array A must
+//>           contain the matrix of coefficients.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>           On entry, LDA specifies the first dimension of A as declared
+//>           in the calling (sub) program. LDA must be at least
+//>           max( 1, m ).
+//> \endverbatim
+//>
+//> \param[in] X
+//> \verbatim
+//>          X is DOUBLE PRECISION array, dimension at least
+//>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
+//>           and at least
+//>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
+//>           Before entry, the incremented array X must contain the
+//>           vector x.
+//> \endverbatim
+//>
+//> \param[in] INCX
+//> \verbatim
+//>          INCX is INTEGER
+//>           On entry, INCX specifies the increment for the elements of
+//>           X. INCX must not be zero.
+//> \endverbatim
+//>
+//> \param[in] BETA
+//> \verbatim
+//>          BETA is DOUBLE PRECISION.
+//>           On entry, BETA specifies the scalar beta. When BETA is
+//>           supplied as zero then Y need not be set on input.
+//> \endverbatim
+//>
+//> \param[in,out] Y
+//> \verbatim
+//>          Y is DOUBLE PRECISION array, dimension at least
+//>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
+//>           and at least
+//>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
+//>           Before entry with BETA non-zero, the incremented array Y
+//>           must contain the vector y. On exit, Y is overwritten by the
+//>           updated vector y.
+//> \endverbatim
+//>
+//> \param[in] INCY
+//> \verbatim
+//>          INCY is INTEGER
+//>           On entry, INCY specifies the increment for the elements of
+//>           Y. INCY must not be zero.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup double_blas_level2
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>  Level 2 Blas routine.
+//>  The vector and matrix arguments are not referenced when N = 0, or M = 0
+//>
+//>  -- Written on 22-October-1986.
+//>     Jack Dongarra, Argonne National Lab.
+//>     Jeremy Du Croz, Nag Central Office.
+//>     Sven Hammarling, Nag Central Office.
+//>     Richard Hanson, Sandia National Labs.
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int dgemv_(char *trans, int *m, int *n, double *alpha,
+	double *a, int *lda, double *x, int *incx, double *beta, double *y,
+	int *incy)
+{
+    // System generated locals
+    int a_dim1, a_offset, i__1, i__2;
+
+    // Local variables
+    int i__, j, ix, iy, jx, jy, kx, ky, info;
+    double temp;
+    int lenx, leny;
+    extern int lsame_(char *, char *);
+    extern /* Subroutine */ int xerbla_(char *, int *);
+
+    //
+    // -- Reference BLAS level2 routine (version 3.7.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //
+    //    Test the input parameters.
+    //
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --x;
+    --y;
+
+    // Function Body
+    info = 0;
+    if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C"))
+	    {
+	info = 1;
+    } else if (*m < 0) {
+	info = 2;
+    } else if (*n < 0) {
+	info = 3;
+    } else if (*lda < max(1,*m)) {
+	info = 6;
+    } else if (*incx == 0) {
+	info = 8;
+    } else if (*incy == 0) {
+	info = 11;
+    }
+    if (info != 0) {
+	xerbla_("DGEMV ", &info);
+	return 0;
+    }
+    //
+    //    Quick return if possible.
+    //
+    if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) {
+	return 0;
+    }
+    //
+    //    Set  LENX  and  LENY, the lengths of the vectors x and y, and set
+    //    up the start points in  X  and  Y.
+    //
+    if (lsame_(trans, "N")) {
+	lenx = *n;
+	leny = *m;
+    } else {
+	lenx = *m;
+	leny = *n;
+    }
+    if (*incx > 0) {
+	kx = 1;
+    } else {
+	kx = 1 - (lenx - 1) * *incx;
+    }
+    if (*incy > 0) {
+	ky = 1;
+    } else {
+	ky = 1 - (leny - 1) * *incy;
+    }
+    //
+    //    Start the operations. In this version the elements of A are
+    //    accessed sequentially with one pass through A.
+    //
+    //    First form  y := beta*y.
+    //
+    if (*beta != 1.) {
+	if (*incy == 1) {
+	    if (*beta == 0.) {
+		i__1 = leny;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = 0.;
+// L10:
+		}
+	    } else {
+		i__1 = leny;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[i__] = *beta * y[i__];
+// L20:
+		}
+	    }
+	} else {
+	    iy = ky;
+	    if (*beta == 0.) {
+		i__1 = leny;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = 0.;
+		    iy += *incy;
+// L30:
+		}
+	    } else {
+		i__1 = leny;
+		for (i__ = 1; i__ <= i__1; ++i__) {
+		    y[iy] = *beta * y[iy];
+		    iy += *incy;
+// L40:
+		}
+	    }
+	}
+    }
+    if (*alpha == 0.) {
+	return 0;
+    }
+    if (lsame_(trans, "N")) {
+	//
+	//       Form  y := alpha*A*x + y.
+	//
+	jx = kx;
+	if (*incy == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = *alpha * x[jx];
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    y[i__] += temp * a[i__ + j * a_dim1];
+// L50:
+		}
+		jx += *incx;
+// L60:
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = *alpha * x[jx];
+		iy = ky;
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    y[iy] += temp * a[i__ + j * a_dim1];
+		    iy += *incy;
+// L70:
+		}
+		jx += *incx;
+// L80:
+	    }
+	}
+    } else {
+	//
+	//       Form  y := alpha*A**T*x + y.
+	//
+	jy = ky;
+	if (*incx == 1) {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = 0.;
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp += a[i__ + j * a_dim1] * x[i__];
+// L90:
+		}
+		y[jy] += *alpha * temp;
+		jy += *incy;
+// L100:
+	    }
+	} else {
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = 0.;
+		ix = kx;
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    temp += a[i__ + j * a_dim1] * x[ix];
+		    ix += *incx;
+// L110:
+		}
+		y[jy] += *alpha * temp;
+		jy += *incy;
+// L120:
+	    }
+	}
+    }
+    return 0;
+    //
+    //    End of DGEMV .
+    //
+} // dgemv_
+
@@ -0,0 +1,186 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DISNAN tests input for NaN.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DISNAN + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/disnan.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/disnan.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/disnan.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      LOGICAL FUNCTION DISNAN( DIN )
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION, INTENT(IN) :: DIN
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DISNAN returns .TRUE. if its argument is NaN, and .FALSE.
+//> otherwise.  To be replaced by the Fortran 2003 intrinsic in the
+//> future.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] DIN
+//> \verbatim
+//>          DIN is DOUBLE PRECISION
+//>          Input to test for NaN.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date June 2017
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+int disnan_(double *din)
+{
+    // System generated locals
+    int ret_val;
+
+    // Local variables
+    extern int dlaisnan_(double *, double *);
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.1) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    June 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    // .. External Functions ..
+    // ..
+    // .. Executable Statements ..
+    ret_val = dlaisnan_(din, din);
+    return ret_val;
+} // disnan_
+
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+//> \brief \b DLAISNAN tests input for NaN by comparing two arguments for inequality.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLAISNAN + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaisnan.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaisnan.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaisnan.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      LOGICAL FUNCTION DLAISNAN( DIN1, DIN2 )
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION, INTENT(IN) :: DIN1, DIN2
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> This routine is not for general use.  It exists solely to avoid
+//> over-optimization in DISNAN.
+//>
+//> DLAISNAN checks for NaNs by comparing its two arguments for
+//> inequality.  NaN is the only floating-point value where NaN != NaN
+//> returns .TRUE.  To check for NaNs, pass the same variable as both
+//> arguments.
+//>
+//> A compiler must assume that the two arguments are
+//> not the same variable, and the test will not be optimized away.
+//> Interprocedural or whole-program optimization may delete this
+//> test.  The ISNAN functions will be replaced by the correct
+//> Fortran 03 intrinsic once the intrinsic is widely available.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] DIN1
+//> \verbatim
+//>          DIN1 is DOUBLE PRECISION
+//> \endverbatim
+//>
+//> \param[in] DIN2
+//> \verbatim
+//>          DIN2 is DOUBLE PRECISION
+//>          Two numbers to compare for inequality.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date June 2017
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+int dlaisnan_(double *din1, double *din2)
+{
+    // System generated locals
+    int ret_val;
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.1) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    June 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    // .. Executable Statements ..
+    ret_val = *din1 != *din2;
+    return ret_val;
+} // dlaisnan_
+
@@ -0,0 +1,184 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DLACPY copies all or part of one two-dimensional array to another.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLACPY + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlacpy.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlacpy.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlacpy.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB )
+//
+//      .. Scalar Arguments ..
+//      CHARACTER          UPLO
+//      INTEGER            LDA, LDB, M, N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLACPY copies all or part of a two-dimensional matrix A to another
+//> matrix B.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] UPLO
+//> \verbatim
+//>          UPLO is CHARACTER*1
+//>          Specifies the part of the matrix A to be copied to B.
+//>          = 'U':      Upper triangular part
+//>          = 'L':      Lower triangular part
+//>          Otherwise:  All of the matrix A
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix A.  M >= 0.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix A.  N >= 0.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension (LDA,N)
+//>          The m by n matrix A.  If UPLO = 'U', only the upper triangle
+//>          or trapezoid is accessed; if UPLO = 'L', only the lower
+//>          triangle or trapezoid is accessed.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>          The leading dimension of the array A.  LDA >= max(1,M).
+//> \endverbatim
+//>
+//> \param[out] B
+//> \verbatim
+//>          B is DOUBLE PRECISION array, dimension (LDB,N)
+//>          On exit, B = A in the locations specified by UPLO.
+//> \endverbatim
+//>
+//> \param[in] LDB
+//> \verbatim
+//>          LDB is INTEGER
+//>          The leading dimension of the array B.  LDB >= max(1,M).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+/* Subroutine */ int dlacpy_(char *uplo, int *m, int *n, double *a, int *lda,
+	double *b, int *ldb)
+{
+    // System generated locals
+    int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
+
+    // Local variables
+    int i__, j;
+    extern int lsame_(char *, char *);
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    b_dim1 = *ldb;
+    b_offset = 1 + b_dim1;
+    b -= b_offset;
+
+    // Function Body
+    if (lsame_(uplo, "U")) {
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__2 = min(j,*m);
+	    for (i__ = 1; i__ <= i__2; ++i__) {
+		b[i__ + j * b_dim1] = a[i__ + j * a_dim1];
+// L10:
+	    }
+// L20:
+	}
+    } else if (lsame_(uplo, "L")) {
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__2 = *m;
+	    for (i__ = j; i__ <= i__2; ++i__) {
+		b[i__ + j * b_dim1] = a[i__ + j * a_dim1];
+// L30:
+	    }
+// L40:
+	}
+    } else {
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__2 = *m;
+	    for (i__ = 1; i__ <= i__2; ++i__) {
+		b[i__ + j * b_dim1] = a[i__ + j * a_dim1];
+// L50:
+	    }
+// L60:
+	}
+    }
+    return 0;
+    //
+    //    End of DLACPY
+    //
+} // dlacpy_
+
@@ -0,0 +1,367 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DCOMBSSQ adds two scaled sum of squares quantities.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DCOMBSSQ( V1, V2 )
+//
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   V1( 2 ), V2( 2 )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2.
+//> That is,
+//>
+//>    V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq
+//>                            + V2_scale**2 * V2_sumsq
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in,out] V1
+//> \verbatim
+//>          V1 is DOUBLE PRECISION array, dimension (2).
+//>          The first scaled sum.
+//>          V1(1) = V1_scale, V1(2) = V1_sumsq.
+//> \endverbatim
+//>
+//> \param[in] V2
+//> \verbatim
+//>          V2 is DOUBLE PRECISION array, dimension (2).
+//>          The second scaled sum.
+//>          V2(1) = V2_scale, V2(2) = V2_sumsq.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date November 2018
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+/* Subroutine */ int dcombssq_(double *v1, double *v2)
+{
+    // System generated locals
+    double d__1;
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    November 2018
+    //
+    //    .. Array Arguments ..
+    //    ..
+    //
+    //=====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    // Parameter adjustments
+    --v2;
+    --v1;
+
+    // Function Body
+    if (v1[1] >= v2[1]) {
+	if (v1[1] != 0.) {
+	    // Computing 2nd power
+	    d__1 = v2[1] / v1[1];
+	    v1[2] += d__1 * d__1 * v2[2];
+	}
+    } else {
+	// Computing 2nd power
+	d__1 = v1[1] / v2[1];
+	v1[2] = v2[2] + d__1 * d__1 * v1[2];
+	v1[1] = v2[1];
+    }
+    return 0;
+    //
+    //    End of DCOMBSSQ
+    //
+} // dcombssq_
+
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+//> \brief \b DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general rectangular matrix.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLANGE + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlange.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlange.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlange.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      DOUBLE PRECISION FUNCTION DLANGE( NORM, M, N, A, LDA, WORK )
+//
+//      .. Scalar Arguments ..
+//      CHARACTER          NORM
+//      INTEGER            LDA, M, N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   A( LDA, * ), WORK( * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLANGE  returns the value of the one norm,  or the Frobenius norm, or
+//> the  infinity norm,  or the  element of  largest absolute value  of a
+//> real matrix A.
+//> \endverbatim
+//>
+//> \return DLANGE
+//> \verbatim
+//>
+//>    DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm'
+//>             (
+//>             ( norm1(A),         NORM = '1', 'O' or 'o'
+//>             (
+//>             ( normI(A),         NORM = 'I' or 'i'
+//>             (
+//>             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
+//>
+//> where  norm1  denotes the  one norm of a matrix (maximum column sum),
+//> normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
+//> normF  denotes the  Frobenius norm of a matrix (square root of sum of
+//> squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] NORM
+//> \verbatim
+//>          NORM is CHARACTER*1
+//>          Specifies the value to be returned in DLANGE as described
+//>          above.
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix A.  M >= 0.  When M = 0,
+//>          DLANGE is set to zero.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix A.  N >= 0.  When N = 0,
+//>          DLANGE is set to zero.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension (LDA,N)
+//>          The m by n matrix A.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>          The leading dimension of the array A.  LDA >= max(M,1).
+//> \endverbatim
+//>
+//> \param[out] WORK
+//> \verbatim
+//>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
+//>          where LWORK >= M when NORM = 'I'; otherwise, WORK is not
+//>          referenced.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup doubleGEauxiliary
+//
+// =====================================================================
+double dlange_(char *norm, int *m, int *n, double *a, int *lda, double *work)
+{
+    // Table of constant values
+    int c__1 = 1;
+
+    // System generated locals
+    int a_dim1, a_offset, i__1, i__2;
+    double ret_val, d__1;
+
+    // Local variables
+    extern /* Subroutine */ int dcombssq_(double *, double *);
+    int i__, j;
+    double sum, ssq[2], temp;
+    extern int lsame_(char *, char *);
+    double value;
+    extern int disnan_(double *);
+    extern /* Subroutine */ int dlassq_(int *, double *, int *, double *,
+	    double *);
+    double colssq[2];
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    //=====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Local Arrays ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    --work;
+
+    // Function Body
+    if (min(*m,*n) == 0) {
+	value = 0.;
+    } else if (lsame_(norm, "M")) {
+	//
+	//       Find max(abs(A(i,j))).
+	//
+	value = 0.;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__2 = *m;
+	    for (i__ = 1; i__ <= i__2; ++i__) {
+		temp = (d__1 = a[i__ + j * a_dim1], abs(d__1));
+		if (value < temp || disnan_(&temp)) {
+		    value = temp;
+		}
+// L10:
+	    }
+// L20:
+	}
+    } else if (lsame_(norm, "O") || *(unsigned char *)norm == '1') {
+	//
+	//       Find norm1(A).
+	//
+	value = 0.;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    sum = 0.;
+	    i__2 = *m;
+	    for (i__ = 1; i__ <= i__2; ++i__) {
+		sum += (d__1 = a[i__ + j * a_dim1], abs(d__1));
+// L30:
+	    }
+	    if (value < sum || disnan_(&sum)) {
+		value = sum;
+	    }
+// L40:
+	}
+    } else if (lsame_(norm, "I")) {
+	//
+	//       Find normI(A).
+	//
+	i__1 = *m;
+	for (i__ = 1; i__ <= i__1; ++i__) {
+	    work[i__] = 0.;
+// L50:
+	}
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__2 = *m;
+	    for (i__ = 1; i__ <= i__2; ++i__) {
+		work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1));
+// L60:
+	    }
+// L70:
+	}
+	value = 0.;
+	i__1 = *m;
+	for (i__ = 1; i__ <= i__1; ++i__) {
+	    temp = work[i__];
+	    if (value < temp || disnan_(&temp)) {
+		value = temp;
+	    }
+// L80:
+	}
+    } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
+	//
+	//       Find normF(A).
+	//       SSQ(1) is scale
+	//       SSQ(2) is sum-of-squares
+	//       For better accuracy, sum each column separately.
+	//
+	ssq[0] = 0.;
+	ssq[1] = 1.;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    colssq[0] = 0.;
+	    colssq[1] = 1.;
+	    dlassq_(m, &a[j * a_dim1 + 1], &c__1, colssq, &colssq[1]);
+	    dcombssq_(ssq, colssq);
+// L90:
+	}
+	value = ssq[0] * sqrt(ssq[1]);
+    }
+    ret_val = value;
+    return ret_val;
+    //
+    //    End of DLANGE
+    //
+} // dlange_
+
@@ -0,0 +1,125 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DLAPY2 returns sqrt(x2+y2).
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLAPY2 + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlapy2.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlapy2.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlapy2.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION   X, Y
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary
+//> overflow.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] X
+//> \verbatim
+//>          X is DOUBLE PRECISION
+//> \endverbatim
+//>
+//> \param[in] Y
+//> \verbatim
+//>          Y is DOUBLE PRECISION
+//>          X and Y specify the values x and y.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date June 2017
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+double dlapy2_(double *x, double *y)
+{
+    // System generated locals
+    double ret_val, d__1;
+
+    // Local variables
+    int x_is_nan__, y_is_nan__;
+    double w, z__, xabs, yabs;
+    extern int disnan_(double *);
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.1) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    June 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    x_is_nan__ = disnan_(x);
+    y_is_nan__ = disnan_(y);
+    if (x_is_nan__) {
+	ret_val = *x;
+    }
+    if (y_is_nan__) {
+	ret_val = *y;
+    }
+    if (! (x_is_nan__ || y_is_nan__)) {
+	xabs = abs(*x);
+	yabs = abs(*y);
+	w = max(xabs,yabs);
+	z__ = min(xabs,yabs);
+	if (z__ == 0.) {
+	    ret_val = w;
+	} else {
+	    // Computing 2nd power
+	    d__1 = z__ / w;
+	    ret_val = w * sqrt(d__1 * d__1 + 1.);
+	}
+    }
+    return ret_val;
+    //
+    //    End of DLAPY2
+    //
+} // dlapy2_
+
@@ -0,0 +1,768 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DGER
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+//
+//      .. Scalar Arguments ..
+//      DOUBLE PRECISION ALPHA
+//      INTEGER INCX,INCY,LDA,M,N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DGER   performs the rank 1 operation
+//>
+//>    A := alpha*x*y**T + A,
+//>
+//> where alpha is a scalar, x is an m element vector, y is an n element
+//> vector and A is an m by n matrix.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>           On entry, M specifies the number of rows of the matrix A.
+//>           M must be at least zero.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>           On entry, N specifies the number of columns of the matrix A.
+//>           N must be at least zero.
+//> \endverbatim
+//>
+//> \param[in] ALPHA
+//> \verbatim
+//>          ALPHA is DOUBLE PRECISION.
+//>           On entry, ALPHA specifies the scalar alpha.
+//> \endverbatim
+//>
+//> \param[in] X
+//> \verbatim
+//>          X is DOUBLE PRECISION array, dimension at least
+//>           ( 1 + ( m - 1 )*abs( INCX ) ).
+//>           Before entry, the incremented array X must contain the m
+//>           element vector x.
+//> \endverbatim
+//>
+//> \param[in] INCX
+//> \verbatim
+//>          INCX is INTEGER
+//>           On entry, INCX specifies the increment for the elements of
+//>           X. INCX must not be zero.
+//> \endverbatim
+//>
+//> \param[in] Y
+//> \verbatim
+//>          Y is DOUBLE PRECISION array, dimension at least
+//>           ( 1 + ( n - 1 )*abs( INCY ) ).
+//>           Before entry, the incremented array Y must contain the n
+//>           element vector y.
+//> \endverbatim
+//>
+//> \param[in] INCY
+//> \verbatim
+//>          INCY is INTEGER
+//>           On entry, INCY specifies the increment for the elements of
+//>           Y. INCY must not be zero.
+//> \endverbatim
+//>
+//> \param[in,out] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension ( LDA, N )
+//>           Before entry, the leading m by n part of the array A must
+//>           contain the matrix of coefficients. On exit, A is
+//>           overwritten by the updated matrix.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>           On entry, LDA specifies the first dimension of A as declared
+//>           in the calling (sub) program. LDA must be at least
+//>           max( 1, m ).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup double_blas_level2
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>  Level 2 Blas routine.
+//>
+//>  -- Written on 22-October-1986.
+//>     Jack Dongarra, Argonne National Lab.
+//>     Jeremy Du Croz, Nag Central Office.
+//>     Sven Hammarling, Nag Central Office.
+//>     Richard Hanson, Sandia National Labs.
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int dger_(int *m, int *n, double *alpha, double *x, int *
+	incx, double *y, int *incy, double *a, int *lda)
+{
+    // System generated locals
+    int a_dim1, a_offset, i__1, i__2;
+
+    // Local variables
+    int i__, j, ix, jy, kx, info;
+    double temp;
+    extern /* Subroutine */ int xerbla_(char *, int *);
+
+    //
+    // -- Reference BLAS level2 routine (version 3.7.0) --
+    // -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //
+    //    Test the input parameters.
+    //
+    // Parameter adjustments
+    --x;
+    --y;
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    // Function Body
+    info = 0;
+    if (*m < 0) {
+	info = 1;
+    } else if (*n < 0) {
+	info = 2;
+    } else if (*incx == 0) {
+	info = 5;
+    } else if (*incy == 0) {
+	info = 7;
+    } else if (*lda < max(1,*m)) {
+	info = 9;
+    }
+    if (info != 0) {
+	xerbla_("DGER  ", &info);
+	return 0;
+    }
+    //
+    //    Quick return if possible.
+    //
+    if (*m == 0 || *n == 0 || *alpha == 0.) {
+	return 0;
+    }
+    //
+    //    Start the operations. In this version the elements of A are
+    //    accessed sequentially with one pass through A.
+    //
+    if (*incy > 0) {
+	jy = 1;
+    } else {
+	jy = 1 - (*n - 1) * *incy;
+    }
+    if (*incx == 1) {
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    if (y[jy] != 0.) {
+		temp = *alpha * y[jy];
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    a[i__ + j * a_dim1] += x[i__] * temp;
+// L10:
+		}
+	    }
+	    jy += *incy;
+// L20:
+	}
+    } else {
+	if (*incx > 0) {
+	    kx = 1;
+	} else {
+	    kx = 1 - (*m - 1) * *incx;
+	}
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    if (y[jy] != 0.) {
+		temp = *alpha * y[jy];
+		ix = kx;
+		i__2 = *m;
+		for (i__ = 1; i__ <= i__2; ++i__) {
+		    a[i__ + j * a_dim1] += x[ix] * temp;
+		    ix += *incx;
+// L30:
+		}
+	    }
+	    jy += *incy;
+// L40:
+	}
+    }
+    return 0;
+    //
+    //    End of DGER  .
+    //
+} // dger_
+
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+//> \brief \b DLARF applies an elementary reflector to a general rectangular matrix.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLARF + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarf.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarf.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarf.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+//
+//      .. Scalar Arguments ..
+//      CHARACTER          SIDE
+//      INTEGER            INCV, LDC, M, N
+//      DOUBLE PRECISION   TAU
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLARF applies a real elementary reflector H to a real m by n matrix
+//> C, from either the left or the right. H is represented in the form
+//>
+//>       H = I - tau * v * v**T
+//>
+//> where tau is a real scalar and v is a real vector.
+//>
+//> If tau = 0, then H is taken to be the unit matrix.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] SIDE
+//> \verbatim
+//>          SIDE is CHARACTER*1
+//>          = 'L': form  H * C
+//>          = 'R': form  C * H
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix C.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix C.
+//> \endverbatim
+//>
+//> \param[in] V
+//> \verbatim
+//>          V is DOUBLE PRECISION array, dimension
+//>                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
+//>                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
+//>          The vector v in the representation of H. V is not used if
+//>          TAU = 0.
+//> \endverbatim
+//>
+//> \param[in] INCV
+//> \verbatim
+//>          INCV is INTEGER
+//>          The increment between elements of v. INCV <> 0.
+//> \endverbatim
+//>
+//> \param[in] TAU
+//> \verbatim
+//>          TAU is DOUBLE PRECISION
+//>          The value tau in the representation of H.
+//> \endverbatim
+//>
+//> \param[in,out] C
+//> \verbatim
+//>          C is DOUBLE PRECISION array, dimension (LDC,N)
+//>          On entry, the m by n matrix C.
+//>          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
+//>          or C * H if SIDE = 'R'.
+//> \endverbatim
+//>
+//> \param[in] LDC
+//> \verbatim
+//>          LDC is INTEGER
+//>          The leading dimension of the array C. LDC >= max(1,M).
+//> \endverbatim
+//>
+//> \param[out] WORK
+//> \verbatim
+//>          WORK is DOUBLE PRECISION array, dimension
+//>                         (N) if SIDE = 'L'
+//>                      or (M) if SIDE = 'R'
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup doubleOTHERauxiliary
+//
+// =====================================================================
+/* Subroutine */ int dlarf_(char *side, int *m, int *n, double *v, int *incv,
+	double *tau, double *c__, int *ldc, double *work)
+{
+    // Table of constant values
+    double c_b4 = 1.;
+    double c_b5 = 0.;
+    int c__1 = 1;
+
+    // System generated locals
+    int c_dim1, c_offset;
+    double d__1;
+
+    // Local variables
+    int i__;
+    int applyleft;
+    extern /* Subroutine */ int dger_(int *, int *, double *, double *, int *,
+	     double *, int *, double *, int *);
+    extern int lsame_(char *, char *);
+    extern /* Subroutine */ int dgemv_(char *, int *, int *, double *, double
+	    *, int *, double *, int *, double *, double *, int *);
+    int lastc, lastv;
+    extern int iladlc_(int *, int *, double *, int *), iladlr_(int *, int *,
+	    double *, int *);
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    // Parameter adjustments
+    --v;
+    c_dim1 = *ldc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+    --work;
+
+    // Function Body
+    applyleft = lsame_(side, "L");
+    lastv = 0;
+    lastc = 0;
+    if (*tau != 0.) {
+	//    Set up variables for scanning V.  LASTV begins pointing to the end
+	//    of V.
+	if (applyleft) {
+	    lastv = *m;
+	} else {
+	    lastv = *n;
+	}
+	if (*incv > 0) {
+	    i__ = (lastv - 1) * *incv + 1;
+	} else {
+	    i__ = 1;
+	}
+	//    Look for the last non-zero row in V.
+	while(lastv > 0 && v[i__] == 0.) {
+	    --lastv;
+	    i__ -= *incv;
+	}
+	if (applyleft) {
+	    //    Scan for the last non-zero column in C(1:lastv,:).
+	    lastc = iladlc_(&lastv, n, &c__[c_offset], ldc);
+	} else {
+	    //    Scan for the last non-zero row in C(:,1:lastv).
+	    lastc = iladlr_(m, &lastv, &c__[c_offset], ldc);
+	}
+    }
+    //    Note that lastc.eq.0 renders the BLAS operations null; no special
+    //    case is needed at this level.
+    if (applyleft) {
+	//
+	//       Form  H * C
+	//
+	if (lastv > 0) {
+	    //
+	    //          w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
+	    //
+	    dgemv_("Transpose", &lastv, &lastc, &c_b4, &c__[c_offset], ldc, &
+		    v[1], incv, &c_b5, &work[1], &c__1);
+	    //
+	    //          C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T
+	    //
+	    d__1 = -(*tau);
+	    dger_(&lastv, &lastc, &d__1, &v[1], incv, &work[1], &c__1, &c__[
+		    c_offset], ldc);
+	}
+    } else {
+	//
+	//       Form  C * H
+	//
+	if (lastv > 0) {
+	    //
+	    //          w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
+	    //
+	    dgemv_("No transpose", &lastc, &lastv, &c_b4, &c__[c_offset], ldc,
+		     &v[1], incv, &c_b5, &work[1], &c__1);
+	    //
+	    //          C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T
+	    //
+	    d__1 = -(*tau);
+	    dger_(&lastc, &lastv, &d__1, &work[1], &c__1, &v[1], incv, &c__[
+		    c_offset], ldc);
+	}
+    }
+    return 0;
+    //
+    //    End of DLARF
+    //
+} // dlarf_
+
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+//> \brief \b ILADLC scans a matrix for its last non-zero column.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download ILADLC + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/iladlc.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/iladlc.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/iladlc.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      INTEGER FUNCTION ILADLC( M, N, A, LDA )
+//
+//      .. Scalar Arguments ..
+//      INTEGER            M, N, LDA
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   A( LDA, * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> ILADLC scans A for its last non-zero column.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix A.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix A.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension (LDA,N)
+//>          The m by n matrix A.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>          The leading dimension of the array A. LDA >= max(1,M).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+int iladlc_(int *m, int *n, double *a, int *lda)
+{
+    // System generated locals
+    int a_dim1, a_offset, ret_val, i__1;
+
+    // Local variables
+    int i__;
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    //    Quick test for the common case where one corner is non-zero.
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    // Function Body
+    if (*n == 0) {
+	ret_val = *n;
+    } else if (a[*n * a_dim1 + 1] != 0. || a[*m + *n * a_dim1] != 0.) {
+	ret_val = *n;
+    } else {
+	//    Now scan each column from the end, returning with the first non-zero.
+	for (ret_val = *n; ret_val >= 1; --ret_val) {
+	    i__1 = *m;
+	    for (i__ = 1; i__ <= i__1; ++i__) {
+		if (a[i__ + ret_val * a_dim1] != 0.) {
+		    return ret_val;
+		}
+	    }
+	}
+    }
+    return ret_val;
+} // iladlc_
+
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+//> \brief \b ILADLR scans a matrix for its last non-zero row.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download ILADLR + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/iladlr.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/iladlr.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/iladlr.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      INTEGER FUNCTION ILADLR( M, N, A, LDA )
+//
+//      .. Scalar Arguments ..
+//      INTEGER            M, N, LDA
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   A( LDA, * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> ILADLR scans A for its last non-zero row.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix A.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix A.
+//> \endverbatim
+//>
+//> \param[in] A
+//> \verbatim
+//>          A is DOUBLE PRECISION array, dimension (LDA,N)
+//>          The m by n matrix A.
+//> \endverbatim
+//>
+//> \param[in] LDA
+//> \verbatim
+//>          LDA is INTEGER
+//>          The leading dimension of the array A. LDA >= max(1,M).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date December 2016
+//
+//> \ingroup OTHERauxiliary
+//
+// =====================================================================
+int iladlr_(int *m, int *n, double *a, int *lda)
+{
+    // System generated locals
+    int a_dim1, a_offset, ret_val, i__1;
+
+    // Local variables
+    int i__, j;
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    December 2016
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    //    Quick test for the common case where one corner is non-zero.
+    // Parameter adjustments
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    // Function Body
+    if (*m == 0) {
+	ret_val = *m;
+    } else if (a[*m + a_dim1] != 0. || a[*m + *n * a_dim1] != 0.) {
+	ret_val = *m;
+    } else {
+	//    Scan up each column tracking the last zero row seen.
+	ret_val = 0;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    i__ = *m;
+	    while(a[max(i__,1) + j * a_dim1] == 0. && i__ >= 1) {
+		--i__;
+	    }
+	    ret_val = max(ret_val,i__);
+	}
+    }
+    return ret_val;
+} // iladlr_
+
@@ -0,0 +1,824 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DLARFB applies a block reflector or its transpose to a general rectangular matrix.
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLARFB + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+//                         T, LDT, C, LDC, WORK, LDWORK )
+//
+//      .. Scalar Arguments ..
+//      CHARACTER          DIRECT, SIDE, STOREV, TRANS
+//      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
+//     $                   WORK( LDWORK, * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLARFB applies a real block reflector H or its transpose H**T to a
+//> real m by n matrix C, from either the left or the right.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] SIDE
+//> \verbatim
+//>          SIDE is CHARACTER*1
+//>          = 'L': apply H or H**T from the Left
+//>          = 'R': apply H or H**T from the Right
+//> \endverbatim
+//>
+//> \param[in] TRANS
+//> \verbatim
+//>          TRANS is CHARACTER*1
+//>          = 'N': apply H (No transpose)
+//>          = 'T': apply H**T (Transpose)
+//> \endverbatim
+//>
+//> \param[in] DIRECT
+//> \verbatim
+//>          DIRECT is CHARACTER*1
+//>          Indicates how H is formed from a product of elementary
+//>          reflectors
+//>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+//>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+//> \endverbatim
+//>
+//> \param[in] STOREV
+//> \verbatim
+//>          STOREV is CHARACTER*1
+//>          Indicates how the vectors which define the elementary
+//>          reflectors are stored:
+//>          = 'C': Columnwise
+//>          = 'R': Rowwise
+//> \endverbatim
+//>
+//> \param[in] M
+//> \verbatim
+//>          M is INTEGER
+//>          The number of rows of the matrix C.
+//> \endverbatim
+//>
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The number of columns of the matrix C.
+//> \endverbatim
+//>
+//> \param[in] K
+//> \verbatim
+//>          K is INTEGER
+//>          The order of the matrix T (= the number of elementary
+//>          reflectors whose product defines the block reflector).
+//>          If SIDE = 'L', M >= K >= 0;
+//>          if SIDE = 'R', N >= K >= 0.
+//> \endverbatim
+//>
+//> \param[in] V
+//> \verbatim
+//>          V is DOUBLE PRECISION array, dimension
+//>                                (LDV,K) if STOREV = 'C'
+//>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+//>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+//>          The matrix V. See Further Details.
+//> \endverbatim
+//>
+//> \param[in] LDV
+//> \verbatim
+//>          LDV is INTEGER
+//>          The leading dimension of the array V.
+//>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+//>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+//>          if STOREV = 'R', LDV >= K.
+//> \endverbatim
+//>
+//> \param[in] T
+//> \verbatim
+//>          T is DOUBLE PRECISION array, dimension (LDT,K)
+//>          The triangular k by k matrix T in the representation of the
+//>          block reflector.
+//> \endverbatim
+//>
+//> \param[in] LDT
+//> \verbatim
+//>          LDT is INTEGER
+//>          The leading dimension of the array T. LDT >= K.
+//> \endverbatim
+//>
+//> \param[in,out] C
+//> \verbatim
+//>          C is DOUBLE PRECISION array, dimension (LDC,N)
+//>          On entry, the m by n matrix C.
+//>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
+//> \endverbatim
+//>
+//> \param[in] LDC
+//> \verbatim
+//>          LDC is INTEGER
+//>          The leading dimension of the array C. LDC >= max(1,M).
+//> \endverbatim
+//>
+//> \param[out] WORK
+//> \verbatim
+//>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
+//> \endverbatim
+//>
+//> \param[in] LDWORK
+//> \verbatim
+//>          LDWORK is INTEGER
+//>          The leading dimension of the array WORK.
+//>          If SIDE = 'L', LDWORK >= max(1,N);
+//>          if SIDE = 'R', LDWORK >= max(1,M).
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date June 2013
+//
+//> \ingroup doubleOTHERauxiliary
+//
+//> \par Further Details:
+// =====================
+//>
+//> \verbatim
+//>
+//>  The shape of the matrix V and the storage of the vectors which define
+//>  the H(i) is best illustrated by the following example with n = 5 and
+//>  k = 3. The elements equal to 1 are not stored; the corresponding
+//>  array elements are modified but restored on exit. The rest of the
+//>  array is not used.
+//>
+//>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+//>
+//>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+//>                   ( v1  1    )                     (     1 v2 v2 v2 )
+//>                   ( v1 v2  1 )                     (        1 v3 v3 )
+//>                   ( v1 v2 v3 )
+//>                   ( v1 v2 v3 )
+//>
+//>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+//>
+//>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+//>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+//>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+//>                   (     1 v3 )
+//>                   (        1 )
+//> \endverbatim
+//>
+// =====================================================================
+/* Subroutine */ int dlarfb_(char *side, char *trans, char *direct, char *
+	storev, int *m, int *n, int *k, double *v, int *ldv, double *t, int *
+	ldt, double *c__, int *ldc, double *work, int *ldwork)
+{
+    // Table of constant values
+    int c__1 = 1;
+    double c_b14 = 1.;
+    double c_b25 = -1.;
+
+    // System generated locals
+    int c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1,
+	    work_offset, i__1, i__2;
+
+    // Local variables
+    int i__, j;
+    extern /* Subroutine */ int dgemm_(char *, char *, int *, int *, int *,
+	    double *, double *, int *, double *, int *, double *, double *,
+	    int *);
+    extern int lsame_(char *, char *);
+    extern /* Subroutine */ int dcopy_(int *, double *, int *, double *, int *
+	    ), dtrmm_(char *, char *, char *, char *, int *, int *, double *,
+	    double *, int *, double *, int *);
+    char transt[1+1]={'\0'};
+
+    //
+    // -- LAPACK auxiliary routine (version 3.7.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    June 2013
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    //    Quick return if possible
+    //
+    // Parameter adjustments
+    v_dim1 = *ldv;
+    v_offset = 1 + v_dim1;
+    v -= v_offset;
+    t_dim1 = *ldt;
+    t_offset = 1 + t_dim1;
+    t -= t_offset;
+    c_dim1 = *ldc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+    work_dim1 = *ldwork;
+    work_offset = 1 + work_dim1;
+    work -= work_offset;
+
+    // Function Body
+    if (*m <= 0 || *n <= 0) {
+	return 0;
+    }
+    if (lsame_(trans, "N")) {
+	*(unsigned char *)transt = 'T';
+    } else {
+	*(unsigned char *)transt = 'N';
+    }
+    if (lsame_(storev, "C")) {
+	if (lsame_(direct, "F")) {
+	    //
+	    //          Let  V =  ( V1 )    (first K rows)
+	    //                    ( V2 )
+	    //          where  V1  is unit lower triangular.
+	    //
+	    if (lsame_(side, "L")) {
+		//
+		//             Form  H * C  or  H**T * C  where  C = ( C1 )
+		//                                                   ( C2 )
+		//
+		//             W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+		//
+		//             W := C1**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
+			     &c__1);
+// L10:
+		}
+		//
+		//             W := W * V1
+		//
+		dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
+			 &v[v_offset], ldv, &work[work_offset], ldwork);
+		if (*m > *k) {
+		    //
+		    //                W := W + C2**T * V2
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
+			    c__[*k + 1 + c_dim1], ldc, &v[*k + 1 + v_dim1],
+			    ldv, &c_b14, &work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T**T  or  W * T
+		//
+		dtrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - V * W**T
+		//
+		if (*m > *k) {
+		    //
+		    //                C2 := C2 - V2 * W**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
+			    v[*k + 1 + v_dim1], ldv, &work[work_offset],
+			    ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
+		}
+		//
+		//             W := W * V1**T
+		//
+		dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
+			v[v_offset], ldv, &work[work_offset], ldwork);
+		//
+		//             C1 := C1 - W**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *n;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
+// L20:
+		    }
+// L30:
+		}
+	    } else if (lsame_(side, "R")) {
+		//
+		//             Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+		//
+		//             W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+		//
+		//             W := C1
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
+			    work_dim1 + 1], &c__1);
+// L40:
+		}
+		//
+		//             W := W * V1
+		//
+		dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
+			 &v[v_offset], ldv, &work[work_offset], ldwork);
+		if (*n > *k) {
+		    //
+		    //                W := W + C2 * V2
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "No transpose", m, k, &i__1, &
+			    c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc, &v[*k +
+			    1 + v_dim1], ldv, &c_b14, &work[work_offset],
+			    ldwork);
+		}
+		//
+		//             W := W * T  or  W * T**T
+		//
+		dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - W * V**T
+		//
+		if (*n > *k) {
+		    //
+		    //                C2 := C2 - W * V2**T
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
+			    work[work_offset], ldwork, &v[*k + 1 + v_dim1],
+			    ldv, &c_b14, &c__[(*k + 1) * c_dim1 + 1], ldc);
+		}
+		//
+		//             W := W * V1**T
+		//
+		dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
+			v[v_offset], ldv, &work[work_offset], ldwork);
+		//
+		//             C1 := C1 - W
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
+// L50:
+		    }
+// L60:
+		}
+	    }
+	} else {
+	    //
+	    //          Let  V =  ( V1 )
+	    //                    ( V2 )    (last K rows)
+	    //          where  V2  is unit upper triangular.
+	    //
+	    if (lsame_(side, "L")) {
+		//
+		//             Form  H * C  or  H**T * C  where  C = ( C1 )
+		//                                                   ( C2 )
+		//
+		//             W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+		//
+		//             W := C2**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
+			    work_dim1 + 1], &c__1);
+// L70:
+		}
+		//
+		//             W := W * V2
+		//
+		dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
+			 &v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
+			ldwork);
+		if (*m > *k) {
+		    //
+		    //                W := W + C1**T * V1
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "No transpose", n, k, &i__1, &c_b14, &
+			    c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
+			    work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T**T  or  W * T
+		//
+		dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - V * W**T
+		//
+		if (*m > *k) {
+		    //
+		    //                C1 := C1 - V1 * W**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("No transpose", "Transpose", &i__1, n, k, &c_b25, &
+			    v[v_offset], ldv, &work[work_offset], ldwork, &
+			    c_b14, &c__[c_offset], ldc);
+		}
+		//
+		//             W := W * V2**T
+		//
+		dtrmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
+			v[*m - *k + 1 + v_dim1], ldv, &work[work_offset],
+			ldwork);
+		//
+		//             C2 := C2 - W**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *n;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
+				work_dim1];
+// L80:
+		    }
+// L90:
+		}
+	    } else if (lsame_(side, "R")) {
+		//
+		//             Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+		//
+		//             W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+		//
+		//             W := C2
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
+			    j * work_dim1 + 1], &c__1);
+// L100:
+		}
+		//
+		//             W := W * V2
+		//
+		dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
+			 &v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
+			ldwork);
+		if (*n > *k) {
+		    //
+		    //                W := W + C1 * V1
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "No transpose", m, k, &i__1, &
+			    c_b14, &c__[c_offset], ldc, &v[v_offset], ldv, &
+			    c_b14, &work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T  or  W * T**T
+		//
+		dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - W * V**T
+		//
+		if (*n > *k) {
+		    //
+		    //                C1 := C1 - W * V1**T
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, &
+			    work[work_offset], ldwork, &v[v_offset], ldv, &
+			    c_b14, &c__[c_offset], ldc);
+		}
+		//
+		//             W := W * V2**T
+		//
+		dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
+			v[*n - *k + 1 + v_dim1], ldv, &work[work_offset],
+			ldwork);
+		//
+		//             C2 := C2 - W
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
+				work_dim1];
+// L110:
+		    }
+// L120:
+		}
+	    }
+	}
+    } else if (lsame_(storev, "R")) {
+	if (lsame_(direct, "F")) {
+	    //
+	    //          Let  V =  ( V1  V2 )    (V1: first K columns)
+	    //          where  V1  is unit upper triangular.
+	    //
+	    if (lsame_(side, "L")) {
+		//
+		//             Form  H * C  or  H**T * C  where  C = ( C1 )
+		//                                                   ( C2 )
+		//
+		//             W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+		//
+		//             W := C1**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1],
+			     &c__1);
+// L130:
+		}
+		//
+		//             W := W * V1**T
+		//
+		dtrmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, &
+			v[v_offset], ldv, &work[work_offset], ldwork);
+		if (*m > *k) {
+		    //
+		    //                W := W + C2**T * V2**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
+			    c__[*k + 1 + c_dim1], ldc, &v[(*k + 1) * v_dim1 +
+			    1], ldv, &c_b14, &work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T**T  or  W * T
+		//
+		dtrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - V**T * W**T
+		//
+		if (*m > *k) {
+		    //
+		    //                C2 := C2 - V2**T * W**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[(
+			    *k + 1) * v_dim1 + 1], ldv, &work[work_offset],
+			    ldwork, &c_b14, &c__[*k + 1 + c_dim1], ldc);
+		}
+		//
+		//             W := W * V1
+		//
+		dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
+			 &v[v_offset], ldv, &work[work_offset], ldwork);
+		//
+		//             C1 := C1 - W**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *n;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
+// L140:
+		    }
+// L150:
+		}
+	    } else if (lsame_(side, "R")) {
+		//
+		//             Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+		//
+		//             W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+		//
+		//             W := C1
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j *
+			    work_dim1 + 1], &c__1);
+// L160:
+		}
+		//
+		//             W := W * V1**T
+		//
+		dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
+			v[v_offset], ldv, &work[work_offset], ldwork);
+		if (*n > *k) {
+		    //
+		    //                W := W + C2 * V2**T
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
+			    c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) *
+			    v_dim1 + 1], ldv, &c_b14, &work[work_offset],
+			    ldwork);
+		}
+		//
+		//             W := W * T  or  W * T**T
+		//
+		dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - W * V
+		//
+		if (*n > *k) {
+		    //
+		    //                C2 := C2 - W * V2
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "No transpose", m, &i__1, k, &
+			    c_b25, &work[work_offset], ldwork, &v[(*k + 1) *
+			    v_dim1 + 1], ldv, &c_b14, &c__[(*k + 1) * c_dim1
+			    + 1], ldc);
+		}
+		//
+		//             W := W * V1
+		//
+		dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
+			 &v[v_offset], ldv, &work[work_offset], ldwork);
+		//
+		//             C1 := C1 - W
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
+// L170:
+		    }
+// L180:
+		}
+	    }
+	} else {
+	    //
+	    //          Let  V =  ( V1  V2 )    (V2: last K columns)
+	    //          where  V2  is unit lower triangular.
+	    //
+	    if (lsame_(side, "L")) {
+		//
+		//             Form  H * C  or  H**T * C  where  C = ( C1 )
+		//                                                   ( C2 )
+		//
+		//             W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+		//
+		//             W := C2**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j *
+			    work_dim1 + 1], &c__1);
+// L190:
+		}
+		//
+		//             W := W * V2**T
+		//
+		dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
+			v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
+			, ldwork);
+		if (*m > *k) {
+		    //
+		    //                W := W + C1**T * V1**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
+			    c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
+			    work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T**T  or  W * T
+		//
+		dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - V**T * W**T
+		//
+		if (*m > *k) {
+		    //
+		    //                C1 := C1 - V1**T * W**T
+		    //
+		    i__1 = *m - *k;
+		    dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
+			    v_offset], ldv, &work[work_offset], ldwork, &
+			    c_b14, &c__[c_offset], ldc);
+		}
+		//
+		//             W := W * V2
+		//
+		dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
+			 &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
+			work_offset], ldwork);
+		//
+		//             C2 := C2 - W**T
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *n;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j *
+				work_dim1];
+// L200:
+		    }
+// L210:
+		}
+	    } else if (lsame_(side, "R")) {
+		//
+		//             Form  C * H  or  C * H'  where  C = ( C1  C2 )
+		//
+		//             W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+		//
+		//             W := C2
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    dcopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
+			    j * work_dim1 + 1], &c__1);
+// L220:
+		}
+		//
+		//             W := W * V2**T
+		//
+		dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
+			v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
+			, ldwork);
+		if (*n > *k) {
+		    //
+		    //                W := W + C1 * V1**T
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
+			    c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
+			    work[work_offset], ldwork);
+		}
+		//
+		//             W := W * T  or  W * T**T
+		//
+		dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
+			t_offset], ldt, &work[work_offset], ldwork);
+		//
+		//             C := C - W * V
+		//
+		if (*n > *k) {
+		    //
+		    //                C1 := C1 - W * V1
+		    //
+		    i__1 = *n - *k;
+		    dgemm_("No transpose", "No transpose", m, &i__1, k, &
+			    c_b25, &work[work_offset], ldwork, &v[v_offset],
+			    ldv, &c_b14, &c__[c_offset], ldc);
+		}
+		//
+		//             W := W * V2
+		//
+		dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
+			 &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
+			work_offset], ldwork);
+		//
+		//             C1 := C1 - W
+		//
+		i__1 = *k;
+		for (j = 1; j <= i__1; ++j) {
+		    i__2 = *m;
+		    for (i__ = 1; i__ <= i__2; ++i__) {
+			c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j *
+				work_dim1];
+// L230:
+		    }
+// L240:
+		}
+	    }
+	}
+    }
+    return 0;
+    //
+    //    End of DLARFB
+    //
+} // dlarfb_
+
@@ -0,0 +1,216 @@
+/*  -- translated by f2c (version 20201020 (for_lapack)). -- */
+
+#include "f2c.h"
+
+//> \brief \b DLARFG generates an elementary reflector (Householder matrix).
+//
+// =========== DOCUMENTATION ===========
+//
+// Online html documentation available at
+//           http://www.netlib.org/lapack/explore-html/
+//
+//> \htmlonly
+//> Download DLARFG + dependencies
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfg.f">
+//> [TGZ]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfg.f">
+//> [ZIP]</a>
+//> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfg.f">
+//> [TXT]</a>
+//> \endhtmlonly
+//
+// Definition:
+// ===========
+//
+//      SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
+//
+//      .. Scalar Arguments ..
+//      INTEGER            INCX, N
+//      DOUBLE PRECISION   ALPHA, TAU
+//      ..
+//      .. Array Arguments ..
+//      DOUBLE PRECISION   X( * )
+//      ..
+//
+//
+//> \par Purpose:
+// =============
+//>
+//> \verbatim
+//>
+//> DLARFG generates a real elementary reflector H of order n, such
+//> that
+//>
+//>       H * ( alpha ) = ( beta ),   H**T * H = I.
+//>           (   x   )   (   0  )
+//>
+//> where alpha and beta are scalars, and x is an (n-1)-element real
+//> vector. H is represented in the form
+//>
+//>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
+//>                     ( v )
+//>
+//> where tau is a real scalar and v is a real (n-1)-element
+//> vector.
+//>
+//> If the elements of x are all zero, then tau = 0 and H is taken to be
+//> the unit matrix.
+//>
+//> Otherwise  1 <= tau <= 2.
+//> \endverbatim
+//
+// Arguments:
+// ==========
+//
+//> \param[in] N
+//> \verbatim
+//>          N is INTEGER
+//>          The order of the elementary reflector.
+//> \endverbatim
+//>
+//> \param[in,out] ALPHA
+//> \verbatim
+//>          ALPHA is DOUBLE PRECISION
+//>          On entry, the value alpha.
+//>          On exit, it is overwritten with the value beta.
+//> \endverbatim
+//>
+//> \param[in,out] X
+//> \verbatim
+//>          X is DOUBLE PRECISION array, dimension
+//>                         (1+(N-2)*abs(INCX))
+//>          On entry, the vector x.
+//>          On exit, it is overwritten with the vector v.
+//> \endverbatim
+//>
+//> \param[in] INCX
+//> \verbatim
+//>          INCX is INTEGER
+//>          The increment between elements of X. INCX > 0.
+//> \endverbatim
+//>
+//> \param[out] TAU
+//> \verbatim
+//>          TAU is DOUBLE PRECISION
+//>          The value tau.
+//> \endverbatim
+//
+// Authors:
+// ========
+//
+//> \author Univ. of Tennessee
+//> \author Univ. of California Berkeley
+//> \author Univ. of Colorado Denver
+//> \author NAG Ltd.
+//
+//> \date November 2017
+//
+//> \ingroup doubleOTHERauxiliary
+//
+// =====================================================================
+/* Subroutine */ int dlarfg_(int *n, double *alpha, double *x, int *incx,
+	double *tau)
+{
+    // System generated locals
+    int i__1;
+    double d__1;
+
+    // Local variables
+    int j, knt;
+    double beta;
+    extern double dnrm2_(int *, double *, int *);
+    extern /* Subroutine */ int dscal_(int *, double *, double *, int *);
+    double xnorm;
+    extern double dlapy2_(double *, double *), dlamch_(char *);
+    double safmin, rsafmn;
+
+    //
+    // -- LAPACK auxiliary routine (version 3.8.0) --
+    // -- LAPACK is a software package provided by Univ. of Tennessee,    --
+    // -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+    //    November 2017
+    //
+    //    .. Scalar Arguments ..
+    //    ..
+    //    .. Array Arguments ..
+    //    ..
+    //
+    // =====================================================================
+    //
+    //    .. Parameters ..
+    //    ..
+    //    .. Local Scalars ..
+    //    ..
+    //    .. External Functions ..
+    //    ..
+    //    .. Intrinsic Functions ..
+    //    ..
+    //    .. External Subroutines ..
+    //    ..
+    //    .. Executable Statements ..
+    //
+    // Parameter adjustments
+    --x;
+
+    // Function Body
+    if (*n <= 1) {
+	*tau = 0.;
+	return 0;
+    }
+    i__1 = *n - 1;
+    xnorm = dnrm2_(&i__1, &x[1], incx);
+    if (xnorm == 0.) {
+	//
+	//       H  =  I
+	//
+	*tau = 0.;
+    } else {
+	//
+	//       general case
+	//
+	d__1 = dlapy2_(alpha, &xnorm);
+	beta = -d_sign(&d__1, alpha);
+	safmin = dlamch_("S") / dlamch_("E");
+	knt = 0;
+	if (abs(beta) < safmin) {
+	    //
+	    //          XNORM, BETA may be inaccurate; scale X and recompute them
+	    //
+	    rsafmn = 1. / safmin;
+L10:
+	    ++knt;
+	    i__1 = *n - 1;
+	    dscal_(&i__1, &rsafmn, &x[1], incx);
+	    beta *= rsafmn;
+	    *alpha *= rsafmn;
+	    if (abs(beta) < safmin && knt < 20) {
+		goto L10;
+	    }
+	    //
+	    //          New BETA is at most 1, at least SAFMIN
+	    //
+	    i__1 = *n - 1;
+	    xnorm = dnrm2_(&i__1, &x[1], incx);
+	    d__1 = dlapy2_(alpha, &xnorm);
+	    beta = -d_sign(&d__1, alpha);
+	}
+	*tau = (beta - *alpha) / beta;
+	i__1 = *n - 1;
+	d__1 = 1. / (*alpha - beta);
+	dscal_(&i__1, &d__1, &x[1], incx);
+	//
+	//       If ALPHA is subnormal, it may lose relative accuracy
+	//
+	i__1 = knt;
+	for (j = 1; j <= i__1; ++j) {
+	    beta *= safmin;
+// L20:
+	}
+	*alpha = beta;
+    }
+    return 0;
+    //
+    //    End of DLARFG
+    //
+} // dlarfg_
+
--- a/Show More
+++ b/Show More