source: branches/dev-cw-2658/RuntimeLibrary/src/jmi/jmi_linear_algebra.c @ 14000

Last change on this file since 14000 was 14000, checked in by Christian Andersson, 5 weeks ago

Minor cleanup. Related to ticket:5868

File size: 4.6 KB
Line 
1/*
2    Copyright (C) 2017 Modelon AB
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License version 3 as published
6    by the Free Software Foundation, or optionally, under the terms of the
7    Common Public License version 1.0 as published by IBM.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License, or the Common Public License, for more details.
13
14    You should have received copies of the GNU General Public License
15    and the Common Public License along with this program.  If not,
16    see <http://www.gnu.org/licenses/> or
17    <http://www.ibm.com/developerworks/library/os-cpl.html/> respectively.
18*/
19
20#include "jmi_linear_algebra.h"
21
22#include <math.h>
23#include <stdlib.h>
24
25/* Computes y = x.*y (element-wise multiplication)*/
26void jmi_linear_algebra_dxemy(jmi_real_t* x, jmi_real_t* y, jmi_int_t N) {
27    jmi_int_t i;
28   
29    for (i = 0; i < N; i++) {
30        y[i] = y[i]*x[i];
31    }
32}
33
34/* Computes a = sqrt( sum( (wi*xi)^2) / N) */
35jmi_real_t jmi_linear_algebra_wrms(jmi_real_t* weights, jmi_real_t* x, jmi_int_t N) {
36    jmi_int_t i;
37    jmi_real_t sum = 0.0, prod;
38   
39    for (i = 0; i < N; i++) {
40        prod = x[i]*weights[i];
41        sum += prod*prod;
42    }
43
44    return sqrt(sum/N);
45}
46
47jmi_real_t jmi_linear_algebra_norm(jmi_real_t* x, jmi_int_t N) {
48    int i = 1;
49   
50    return dnrm2_(&N, x, &i);
51}
52
53jmi_real_t jmi_linear_algebra_dgecon(jmi_real_t* LU, jmi_real_t A_norm, jmi_int_t N, char norm_type) {
54    jmi_real_t rcond = 0.0;
55    int info = 0;
56    jmi_real_t *work = (jmi_real_t*)calloc(4*N, sizeof(jmi_real_t));
57    jmi_int_t *iwork = (jmi_int_t*)calloc(4*N, sizeof(jmi_int_t));
58   
59    dgecon_(&norm_type, &N, LU, &N, &A_norm, &rcond, work, iwork, &info);
60   
61    free(work);
62    free(iwork);
63   
64    return rcond;
65}
66
67/* Computes matrix norms */
68jmi_real_t jmi_linear_algebra_dlange(jmi_real_t* A, jmi_int_t N, char norm_type) {
69    jmi_real_t *work = 0;
70    jmi_real_t norm;
71   
72    if (norm_type == 'I' || norm_type == 'i') { work = (jmi_real_t*)calloc(N, sizeof(jmi_real_t)); }
73   
74    norm = dlange_(&norm_type, &N, &N, A, &N, work); 
75   
76    if (norm_type == 'I' || norm_type == 'i') { free(work); }
77   
78    return norm;
79}
80
81/* Computes y = ax + y */
82void jmi_linear_algebra_daxpy(jmi_real_t a, jmi_real_t* x, jmi_real_t* y, jmi_int_t N) {
83    int i = 1;
84   
85    daxpy_(&N, &a, x, &i, y, &i);
86}
87
88/* Computes a = x^T y */
89jmi_real_t jmi_linear_algebra_ddot(jmi_real_t* x, jmi_real_t* y, jmi_int_t N) {
90    int i = 1;
91   
92    return ddot_(&N, x, &i, y, &i);
93}
94
95/* Computes  y = alpha*A*x + beta*y */
96void jmi_linear_algebra_dgemv(jmi_real_t a, jmi_real_t* A, jmi_real_t* x, jmi_real_t b, jmi_real_t* y, jmi_int_t N, jmi_int_t trans) {
97    char trans_char = trans? 'T':'N'; /* No transposition */
98    int i = 1;
99   
100    dgemv_(&trans_char, &N, &N, &a, A, &N, x, &i, &b, y, &i);
101}
102
103/* Computes A = alpha*x*y**T + A */
104void jmi_linear_algebra_dger(jmi_real_t a, jmi_real_t* x, jmi_real_t *y, jmi_real_t *A, jmi_int_t N) {
105    int i = 1;
106   
107    dger_(&N, &N, &a, x, &i, y, &i, A, &N);
108}
109
110/* Compytes z = ax + by */
111void jmi_linear_algebra_daxpby(jmi_real_t a, jmi_real_t* x, jmi_real_t b, jmi_real_t* y, jmi_real_t* z, jmi_int_t N) {
112    int i;
113   
114    if (a==b) {
115        for (i = 0; i < N; i++)
116            z[i] = a*(x[i]+y[i]);
117    } else if (a==-b) {
118        for (i = 0; i < N; i++)
119            z[i] = a*(x[i]-y[i]);
120    } else {
121        for (i = 0; i < N; i++)
122            z[i] = a*x[i]+b*y[i];
123    }
124}
125
126/* Find the index of the max absolute value */
127jmi_int_t jmi_linear_algebra_idamax(jmi_real_t *x, jmi_int_t N) {
128    int i = 0, j = 0;
129
130    jmi_real_t cmax = JMI_ABS(x[i]);
131    jmi_real_t tmp;
132    for(i=1; i < N; i++) {
133        tmp = JMI_ABS(x[i]);
134        if(tmp > cmax) {
135            cmax = tmp;
136            j=i;
137        }
138    }
139    return j;
140    /*return idamax_(&N, x, &i) - 1; */ /* Compensate for Fortran indexing */
141}
142
143/* Perform LU factorization using Lapack */
144jmi_int_t jmi_linear_algebra_LU_factorize(jmi_real_t* A, jmi_int_t* pivots, jmi_int_t N) {
145    int info = 0;
146
147    dgetrf_(&N, &N, A, &N, pivots, &info);
148
149    return info;
150}
151
152/* Solve with an LU factorized matrix using Lapack */
153jmi_int_t jmi_linear_algebra_LU_solve(jmi_real_t* LU, jmi_int_t* pivots, jmi_real_t* x, jmi_int_t N) {
154    int info = 0, i = 1;
155    char trans = 'N'; /* No transposition */
156   
157    /* Back-solve and get solution in x */
158    dgetrs_(&trans, &N, &i, LU, &N, pivots, x, &N, &info);
159
160    return info;
161}
Note: See TracBrowser for help on using the repository browser.