aoclsparse_?trsm() - 5.2 English - 68552

  1/* ************************************************************************
  2 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <iomanip>
 27#include <iostream>
 28#include <limits>
 29#include <vector>
 30
 31int main(void)
 32{
 33    std::cout << "------------------------------------------------" << std::endl
 34              << " Triangle Solve with Multiple Right Hand Sides" << std::endl
 35              << " sample program for real data type" << std::endl
 36              << "------------------------------------------------" << std::endl
 37              << std::endl;
 38
 39    /*
 40     * This example illustrates how to solve two triangular
 41     * linear system of equations. The real matrix used is
 42     *     | 1  2  0  0 |
 43     * A = | 3  1  2  0 |
 44     *     | 0  3  1  2 |
 45     *     | 0  0  3  1 |
 46     *
 47     * The first linear system is L X = alpha * B with L = tril(A), X and B
 48     * are two dense matrices of size 4x2 stored in column-major format.
 49     *
 50     * The second linear system is U X = alpha * B with U = triu(A), X and B
 51     * are two dense matrices of size 4x2 stored in row-major format.
 52     *
 53     */
 54
 55    // Create a tri-diagonal matrix in CSR format
 56    const aoclsparse_int n = 4, m = 4, nnz = 10;
 57    aoclsparse_int       icrow[5] = {0, 2, 5, 8, 10};
 58    aoclsparse_int       icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
 59    double               aval[10] = {1, 2, 3, 1, 2, 3, 1, 2, 3, 1};
 60    aoclsparse_status    status;
 61    const aoclsparse_int k = 2;
 62    aoclsparse_int       ldb, ldx;
 63    bool                 oki, ok;
 64    double               tol = 20.0 * std::numeric_limits<double>::epsilon();
 65
 66    // Create aoclsparse matrix and its descriptor
 67    aoclsparse_matrix     A;
 68    aoclsparse_index_base base = aoclsparse_index_base_zero;
 69    aoclsparse_order      order;
 70    aoclsparse_mat_descr  descr_a;
 71    aoclsparse_operation  trans;
 72    status = aoclsparse_create_dcsr(&A, base, m, n, nnz, icrow, icol, aval);
 73    if(status != aoclsparse_status_success)
 74    {
 75        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
 76                  << std::endl;
 77        return 3;
 78    }
 79    aoclsparse_create_mat_descr(&descr_a);
 80
 81    /* Case 1: Solving the lower triangular system L^T X = alpha*B.
 82     *     |1, -1|
 83     * B = |2, -2| stored in column-major layout
 84     *     |3, -6| k = 2 and ldb = 4
 85     *     |4, -8|
 86     *
 87     * Linear system L X = B, with solution matrix
 88     *     | 86,-167|
 89     * X = |-29,  56| stored in column-major layout
 90     *     |  9, -18| k = 2 and ldx = 4
 91     *     | -4,   8|
 92     */
 93
 94    double alpha = -1.0;
 95    double X[m * k];
 96
 97    double              B[m * k] = {1, 2, 3, 4, -1, -2, -6, -8};
 98    std::vector<double> XRef;
 99    XRef.assign({86, -29, 9, -4, -167, 56, -18, 8});
100
101    // Indicate to only access the lower triangular part of matrix A
102    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
103    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
104    trans = aoclsparse_operation_transpose;
105
106    order  = aoclsparse_order_column;
107    ldb    = m;
108    ldx    = m;
109    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
110    if(status != aoclsparse_status_success)
111    {
112        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
113                  << std::endl;
114        return 3;
115    }
116    status = aoclsparse_optimize(A);
117    if(status != aoclsparse_status_success)
118    {
119        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
120                  << std::endl;
121        return 3;
122    }
123    // Solve
124    status = aoclsparse_dtrsm(trans, alpha, A, descr_a, order, &B[0], k, ldb, &X[0], ldx);
125    if(status != aoclsparse_status_success)
126    {
127        std::cerr << "Error returned from aoclsparse_dtrsm, status = " << status << "."
128                  << std::endl;
129        return 3;
130    }
131
132    // Print and check the result
133    std::cout << "Solving L^T X = alpha B, where  L=tril(A), and" << std::endl;
134    std::cout << "  X and B are dense rectangular martices (column-major layout)" << std::endl;
135    std::cout << "  Solution matrix X = " << std::endl;
136    std::cout << std::fixed;
137    std::cout.precision(1);
138    ok = true;
139    size_t idx;
140    for(int row = 0; row < m; ++row)
141    {
142        for(int col = 0; col < k; ++col)
143        {
144            idx = row + col * ldx;
145            oki = std::abs(X[idx] - XRef[idx]) <= tol;
146            std::cout << X[idx] << (oki ? "  " : "! ");
147            ok &= oki;
148        }
149        std::cout << std::endl;
150    }
151    std::cout << std::endl;
152
153    /* Case 2: Solving the lower triangular system U X = alpha*B.
154     * In this case we "reinterpret" B as a col-major layout, having the
155     * form
156     *     | 1,  2|
157     * B = | 3,  4| stored in row-major layout
158     *     |-1, -2| k = 2 and ldb = 2
159     *     |-6, -8|
160     *
161     * Linear system L X = B, with solution matrix
162     *     | 39,  50|
163     * X = |-19, -24| stored in row-major layout
164     *     | 11,  14| k = 2 and ldx = 2
165     *     | -6,  -8|
166     */
167
168    // Store same XRef in row-major format
169    XRef.assign({39, 50, -19, -24, 11, 14, -6, -8});
170    alpha = 1.0;
171    ldb   = k;
172    ldx   = k;
173    // Indicate to use only the conjugate transpose of the upper part of A
174    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
175    trans  = aoclsparse_operation_none;
176    order  = aoclsparse_order_row;
177    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
178    if(status != aoclsparse_status_success)
179    {
180        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
181                  << std::endl;
182        return 3;
183    }
184    status = aoclsparse_optimize(A);
185    if(status != aoclsparse_status_success)
186    {
187        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
188                  << std::endl;
189        return 3;
190    }
191    // Solve
192    status = aoclsparse_dtrsm(trans, alpha, A, descr_a, order, &B[0], k, ldb, &X[0], ldx);
193    if(status != aoclsparse_status_success)
194    {
195        std::cerr << "Error returned from aoclsparse_dtrsm, status = " << status << "."
196                  << std::endl;
197        return 3;
198    }
199    // Print and check the result
200    std::cout << "Solving U X = alpha B, where  U=triu(A), and" << std::endl;
201    std::cout << "  X and B are dense rectangular martices (row-major layout)" << std::endl;
202    std::cout << "  Solution matrix X = " << std::endl;
203    for(int row = 0; row < m; ++row)
204    {
205        for(int col = 0; col < k; col++)
206        {
207            idx = row * ldx + col;
208            oki = std::abs(X[idx] - XRef[idx]) <= tol;
209            std::cout << X[idx] << (oki ? " " : "!  ");
210            ok &= oki;
211        }
212        std::cout << std::endl;
213    }
214    std::cout << std::endl;
215
216    // Destroy the aoclsparse memory
217    aoclsparse_destroy_mat_descr(descr_a);
218    aoclsparse_destroy(&A);
219
220    return ok ? 0 : 5;
221}

  1/* ************************************************************************
  2 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <complex>
 27#include <iomanip>
 28#include <iostream>
 29#include <limits>
 30#include <vector>
 31
 32int main(void)
 33{
 34    std::cout << "------------------------------------------------" << std::endl
 35              << " Triangle Solve with Multiple Right Hand Sides" << std::endl
 36              << " sample program for complex data type" << std::endl
 37              << "------------------------------------------------" << std::endl
 38              << std::endl;
 39
 40    /*
 41     * This example illustrates how to solve two triangular
 42     * linear system of equations. The complex matrix used is
 43     *     | 1+3i  2+5i     0     0 |
 44     * A = | 3     1+2i  2-2i     0 |
 45     *     |    0     3     1  2+3i |
 46     *     |    0     0  3+1i  1+1i |
 47     *
 48     * The first linear system is L X = alpha * B with L = tril(A), X and B
 49     * are two dense matrices of size 4x2 stored in column-major format.
 50     *
 51     * The second linear system is U^H X = alpha * B with U = triu(A), X and B
 52     * are two dense matrices of size 4x2 stored in row-major format.
 53     *
 54     */
 55
 56    // Create a tri-diagonal matrix in CSR format
 57    const aoclsparse_int                   n = 4, m = 4, nnz = 10;
 58    aoclsparse_int                         icrow[5] = {0, 2, 5, 8, 10};
 59    aoclsparse_int                         icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
 60    std::vector<aoclsparse_double_complex> aval;
 61    // clang-format off
 62    aval.assign({{1, 3}, {2, 5}, {3, 0}, {1, 2}, {2, -2}, {3, 0},
 63                 {1, 0}, {2, 3}, {3, 1}, {1, 1}});
 64    // clang-format on
 65    aoclsparse_status status;
 66    aoclsparse_int    ldb, ldx;
 67    bool              oki, ok;
 68    double            tol = 20 * std::numeric_limits<double>::epsilon();
 69
 70    aoclsparse_matrix     A;
 71    const aoclsparse_int  k     = 2; // number of columns of X and B
 72    aoclsparse_index_base base  = aoclsparse_index_base_zero;
 73    aoclsparse_order      order = aoclsparse_order_column;
 74    aoclsparse_mat_descr  descr_a;
 75    aoclsparse_operation  trans = aoclsparse_operation_none;
 76    status = aoclsparse_create_zcsr(&A, base, m, n, nnz, icrow, icol, aval.data());
 77    if(status != aoclsparse_status_success)
 78    {
 79        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
 80                  << std::endl;
 81        return 3;
 82    }
 83    aoclsparse_create_mat_descr(&descr_a);
 84
 85    /* Case 1: Solving the lower triangular system L X = B.
 86     *     |-2+4i,   -3+i|
 87     * B = |3+13i,  4+11i| stored in column-major layout
 88     *     |16+7i,  16+2i| k = 2 and ldb = 4
 89     *     |15+11i,10+16i|
 90     *
 91     * Linear system L X = B, with solution matrix
 92     *     |1+i,     i|
 93     * X = |4+2i,    4| stored in column-major layout
 94     *     |4+i,  4+2i| k = 2 and ldx = 4
 95     *     |4,    3+3i|
 96     */
 97
 98    aoclsparse_double_complex alpha = {1.0, 0.};
 99    aoclsparse_double_complex X[m * k];
100
101    std::vector<aoclsparse_double_complex> B;
102    B.assign({{-2, 4}, {3, 13}, {16, 7}, {15, 11}, {-3, 1}, {4, 11}, {16, 2}, {10, 16}});
103    std::vector<aoclsparse_double_complex> XRef;
104    XRef.assign({{1, 1}, {4, 2}, {4, 1}, {4, 0}, {0, 1}, {4, 0}, {4, 2}, {3, 3}});
105
106    // indicate to only access the lower triangular part of matrix A
107    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
108    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
109
110    // Prepare and call solver
111    order  = aoclsparse_order_column;
112    ldb    = m;
113    ldx    = m;
114    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
115    if(status != aoclsparse_status_success)
116    {
117        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
118                  << std::endl;
119        return 3;
120    }
121    status = aoclsparse_optimize(A);
122    if(status != aoclsparse_status_success)
123    {
124        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
125                  << std::endl;
126        return 3;
127    }
128
129    // Solve
130    status = aoclsparse_ztrsm(trans, alpha, A, descr_a, order, B.data(), k, ldb, &X[0], ldx);
131    if(status != aoclsparse_status_success)
132    {
133        std::cerr << "Error from aoclsparse_ztrsm, status = " << status << std::endl;
134        return 3;
135    }
136
137    // Print and check the result
138    std::cout << "Solving L X = alpha B, where  L=tril(A), and" << std::endl;
139    std::cout << "  X and B are dense rectangular martices (column-major layout)" << std::endl;
140    std::cout << "  Solution matrix X = " << std::endl;
141    std::cout << std::fixed;
142    std::cout.precision(1);
143    ok = true;
144    size_t idx;
145    for(int row = 0; row < m; ++row)
146    {
147        for(int col = 0; col < k; ++col)
148        {
149            idx = row + col * ldx;
150            oki = std::abs(X[idx].real - XRef[idx].real) <= tol
151                  && std::abs(X[idx].imag - XRef[idx].imag) <= tol;
152            std::cout << "(" << X[idx].real << ", " << X[idx].imag << "i) " << (oki ? "  " : "! ");
153            ok &= oki;
154        }
155        std::cout << std::endl;
156    }
157    std::cout << std::endl;
158
159    /* Case 2: Solving the lower triangular system U^H X = alpha*B.
160     *     |2+4i,    -1+3i|
161     * B = |9+15i,    6+9i| stored in row-major layout
162     *     |-13+8i,-10+12i| k = 2 and ldx = 2
163     *     |14+15i,  8+20i|
164     *
165     * Linear system L X = B, with solution matrix
166     *     |1+i,     i|
167     * X = |4+2i,    4| stored in row-major layout
168     *     |4+i,  4+2i| k = 2 and ldx = 2
169     *     |4,    3+3i|
170     */
171
172    B.assign({{2, 4}, {-1, 3}, {9, 15}, {6, 9}, {-13, 8}, {-10, 12}, {14, 15}, {8, 20}});
173    // Store same XRef in row-major format
174    XRef.assign({{1, 1}, {0, 1}, {4, 2}, {4, 0}, {4, 1}, {4, 2}, {4, 0}, {3, 3}});
175    alpha = {0, -1};
176    ldb   = k;
177    ldx   = k;
178    // Indicate to use only the conjugate transpose of the upper part of A
179    trans = aoclsparse_operation_conjugate_transpose;
180    order = aoclsparse_order_row;
181    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
182    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
183    if(status != aoclsparse_status_success)
184    {
185        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
186                  << std::endl;
187        return 3;
188    }
189    status = aoclsparse_optimize(A);
190    if(status != aoclsparse_status_success)
191    {
192        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
193                  << std::endl;
194        return 3;
195    }
196    // Solve
197    status = aoclsparse_ztrsm(trans, alpha, A, descr_a, order, B.data(), k, ldb, &X[0], ldx);
198    if(status != aoclsparse_status_success)
199    {
200        std::cerr << "Error from aoclsparse_ztrsm, status = " << status << std::endl;
201        return 3;
202    }
203    // Print and check the result
204    std::cout << "Solving U^H X = alpha B, where  U=triu(A), and" << std::endl;
205    std::cout << "  X and B are dense rectangular martices (row-major layout)" << std::endl;
206    std::cout << "  Solution matrix X = " << std::endl;
207    for(int row = 0; row < m; ++row)
208    {
209        for(int col = 0; col < k; col++)
210        {
211            idx = row * ldx + col;
212            oki = std::abs(X[idx].real - XRef[idx].real) <= tol
213                  && std::abs(X[idx].imag - XRef[idx].imag) <= tol;
214            std::cout << "(" << X[idx].real << ", " << X[idx].imag << "i) " << (oki ? "  " : "! ");
215            ok &= oki;
216        }
217        std::cout << std::endl;
218    }
219    std::cout << std::endl;
220    // Destroy the aoclsparse memory
221    aoclsparse_destroy_mat_descr(descr_a);
222    aoclsparse_destroy(&A);
223
224    return ok ? 0 : 5;
225}