1 /* Siconos is a program dedicated to modeling, simulation and control
2 * of non smooth dynamical systems.
3 *
4 * Copyright 2021 INRIA.
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License");
7 * you may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 */
18
19
20
21 #include <assert.h> // for assert
22 #include <stdio.h> // for fprintf, stderr
23 #include <stdlib.h> // for exit, EXIT_FAILURE
24 #include "CSparseMatrix_internal.h" // for CSparseMatrix, CS_INT
25 #include "NM_conversions.h"
26 #include "SiconosConfig.h" // for WITH_MKL_SPBLAS // IWYU pragma: keep
27
28 #ifdef WITH_MKL_SPBLAS
29 #include "tlsdef.h"
30 #include "MKL_common.h"
31 typedef void (*mkl_dcsrcoo_t)(const __INT_T *job, const __INT_T *n, double *acsr, __INT_T *ja, __INT_T *ia, __INT_T *nnz, double *acoo, __INT_T *rowind, __INT_T *colind, __INT_T *info);
32 typedef void (*mkl_dcsrcsc_t)(const __INT_T *job, const __INT_T *n, double *acsr, __INT_T *ja, __INT_T *ia, double *acsc, __INT_T *ja1, __INT_T *ia1, __INT_T *info);
33
34 tlsvar mkl_dcsrcoo_t mkl_dcsrcoo_p = NULL;
35 tlsvar mkl_dcsrcsc_t mkl_dcsrcsc_p = NULL;
36
37 #endif
38
NM_csc_to_triplet(CSparseMatrix * csc)39 CSparseMatrix* NM_csc_to_triplet(CSparseMatrix* csc)
40 {
41 assert(csc);
42 CSparseMatrix* triplet = cs_spalloc(csc->m, csc->n, csc->p[csc->n], 1, 1);
43
44 CS_INT* Ap = csc->p;
45 CS_INT* Ai = csc->i;
46 double* val = csc->x;
47 for(CS_INT j = 0; j < csc->n; ++j)
48 {
49 for(CS_INT i = Ap[j]; i < Ap[j+1]; ++i)
50 {
51 CSparseMatrix_entry(triplet, Ai[i], j, val[i]);
52 }
53 }
54 return triplet;
55 }
56
NM_csc_to_half_triplet(CSparseMatrix * csc)57 CSparseMatrix* NM_csc_to_half_triplet(CSparseMatrix* csc)
58 {
59 assert(csc);
60 CSparseMatrix* triplet = cs_spalloc(csc->m, csc->n, csc->p[csc->n], 1, 1);
61 if(!triplet) return (cs_done(triplet, NULL, NULL, 0)) ;
62 CS_INT* Ap = csc->p;
63 CS_INT* Ai = csc->i;
64 double* val = csc->x;
65 for(CS_INT j = 0; j < csc->n; ++j)
66 {
67 for(CS_INT i = Ap[j]; i < Ap[j+1]; ++i)
68 {
69 CSparseMatrix_symmetric_entry(triplet, Ai[i], j, val[i]);
70 }
71 }
72 return triplet;
73 }
74
NM_triplet_to_csr(CSparseMatrix * triplet)75 CSparseMatrix* NM_triplet_to_csr(CSparseMatrix* triplet)
76 {
77 #ifdef WITH_MKL_SPBLAS
78 assert(triplet);
79 CHECK_MKL(load_mkl_function("mkl_dcsrcoo", (void**)&mkl_dcsrcoo_p));
80 if(triplet->m != triplet->n)
81 {
82 fprintf(stderr, "NM_triplet_to_csr :: the matrix has to be square\n");
83 exit(EXIT_FAILURE);
84 }
85 CSparseMatrix* csr = cs_spalloc(NSM_NROW_CSR(triplet), NSM_NCOL_CSR(triplet), triplet->nz, 1, 0);
86 assert(csr);
87 csr->nz = -2;
88
89 CS_INT n = csr->n;
90 CS_INT job[6] = {0};
91 CS_INT info = 0;
92 job[0] = 2;
93 (*mkl_dcsrcoo_p)(job, &n, csr->x, csr->i, csr->p, &(triplet->nz), triplet->x, triplet->i, triplet->p, &info);
94
95 return csr;
96 #else
97 fprintf(stderr, "NM_triplet_to_csr :: MKL not enabled\n");
98 exit(EXIT_FAILURE);
99 #endif
100 }
101
NM_csr_to_triplet(CSparseMatrix * csr)102 CSparseMatrix* NM_csr_to_triplet(CSparseMatrix* csr)
103 {
104 #ifdef WITH_MKL_SPBLAS
105 assert(csr);
106 CHECK_MKL(load_mkl_function("mkl_dcsrcoo", (void**)&mkl_dcsrcoo_p));
107 if(csr->m != csr->n)
108 {
109 fprintf(stderr, "NM_csr_to_triplet :: the matrix has to be square\n");
110 exit(EXIT_FAILURE);
111 }
112 CSparseMatrix* triplet = cs_spalloc(csr->m, csr->n, csr->p[csr->m], 1, 1);
113 assert(triplet);
114
115 CS_INT n = csr->n;
116 CS_INT job[6] = {0};
117 job[4] = csr->p[csr->m];
118 job[5] = 3;
119 CS_INT info = 0;
120 (*mkl_dcsrcoo_p)(job, &n, csr->x, csr->i, csr->p, &(csr->p[csr->m]), triplet->x, triplet->i, triplet->p, &info);
121 triplet->nz = csr->p[csr->m];
122
123 return triplet;
124 #else
125 fprintf(stderr, "NM_csr_to_triplet :: MKL not enabled\n");
126 exit(EXIT_FAILURE);
127 #endif
128 }
129
NM_csc_to_csr(CSparseMatrix * csc)130 CSparseMatrix* NM_csc_to_csr(CSparseMatrix* csc)
131 {
132 #ifdef WITH_MKL_SPBLAS
133 assert(csc);
134 CHECK_MKL(load_mkl_function("mkl_dcsrcsc", (void**)&mkl_dcsrcsc_p));
135 if(csc->m != csc->n)
136 {
137 fprintf(stderr, "NM_csc_to_csr :: the matrix has to be square\n");
138 exit(EXIT_FAILURE);
139 }
140 CSparseMatrix* csr = cs_spalloc(NSM_NROW_CSR(csc), NSM_NCOL_CSR(csc), csc->p[csc->n], 1, 0);
141 assert(csr);
142 csr->nz = -2;
143
144 CS_INT n = csr->n;
145 CS_INT job[6] = {0};
146 CS_INT info = 0;
147 job[0] = 1;
148 job[5] = 1;
149 (*mkl_dcsrcsc_p)(job, &n, csr->x, csr->i, csr->p, csc->x, csc->i, csc->p, &info);
150
151 return csr;
152 #else
153 fprintf(stderr, "NM_csc_to_csr :: MKL not enabled\n");
154 exit(EXIT_FAILURE);
155 #endif
156 }
157
NM_csr_to_csc(CSparseMatrix * csr)158 CSparseMatrix* NM_csr_to_csc(CSparseMatrix* csr)
159 {
160 #ifdef WITH_MKL_SPBLAS
161 assert(csr);
162 CHECK_MKL(load_mkl_function("mkl_dcsrcsc", (void**)&mkl_dcsrcsc_p));
163 if(csr->m != csr->n)
164 {
165 fprintf(stderr, "NM_csr_to_csc :: the matrix has to be square\n");
166 exit(EXIT_FAILURE);
167 }
168 assert(csr);
169 CSparseMatrix* csc = cs_spalloc(csr->m, csr->n, csr->nzmax, 1, 0);
170 assert(csc);
171
172 CS_INT n = csr->n;
173 CS_INT job[6] = {0};
174 job[5] = 1;
175 CS_INT info = 0;
176 (*mkl_dcsrcsc_p)(job, &n, csr->x, csr->i, csr->p, csc->x, csc->i, csc->p, &info);
177
178 return csc;
179 #else
180 fprintf(stderr, "NM_csr_to_csc :: MKL not enabled\n");
181 exit(EXIT_FAILURE);
182 #endif
183 }
184