source: trunk/RuntimeLibrary/src/jmi/jmi_linear_algebra.c @ 9915

Last change on this file since 9915 was 9915, checked in by aramle, 3 years ago

#5370 Fixed typo.

File size: 4.3 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 a = sqrt( sum( (wi*xi)^2) / N) */
26jmi_real_t jmi_linear_algebra_wrms(jmi_real_t* weights, jmi_real_t* x, jmi_int_t N) {
27    jmi_int_t i;
28    jmi_real_t sum = 0.0, prod;
29   
30    for (i = 0; i < N; i++) {
31        prod = x[i]*weights[i];
32        sum += prod*prod;
33    }
34
35    return sqrt(sum/N);
36}
37
38jmi_real_t jmi_linear_algebra_norm(jmi_real_t* x, jmi_int_t N) {
39    int i = 1;
40   
41    return dnrm2_(&N, x, &i);
42}
43
44jmi_real_t jmi_linear_algebra_dgecon(jmi_real_t* LU, jmi_real_t A_norm, jmi_int_t N, char norm_type) {
45    jmi_real_t rcond = 0.0;
46    int info = 0;
47    jmi_real_t *work = (jmi_real_t*)calloc(4*N, sizeof(jmi_real_t));
48    jmi_int_t *iwork = (jmi_int_t*)calloc(4*N, sizeof(jmi_int_t));
49   
50    dgecon_(&norm_type, &N, LU, &N, &A_norm, &rcond, work, iwork, &info);
51   
52    free(work);
53    free(iwork);
54   
55    return rcond;
56}
57
58/* Computes matrix norms */
59jmi_real_t jmi_linear_algebra_dlange(jmi_real_t* A, jmi_int_t N, char norm_type) {
60    jmi_real_t *work = 0;
61    jmi_real_t norm;
62   
63    if (norm_type == 'I' || norm_type == 'i') { work = (jmi_real_t*)calloc(N, sizeof(jmi_real_t)); }
64   
65    norm = dlange_(&norm_type, &N, &N, A, &N, work); 
66   
67    if (norm_type == 'I' || norm_type == 'i') { free(work); }
68   
69    return norm;
70}
71
72/* Computes y = ax + y */
73void jmi_linear_algebra_daxpy(jmi_real_t a, jmi_real_t* x, jmi_real_t* y, jmi_int_t N) {
74    int i = 1;
75   
76    daxpy_(&N, &a, x, &i, y, &i);
77}
78
79/* Computes a = x^T y */
80jmi_real_t jmi_linear_algebra_ddot(jmi_real_t* x, jmi_real_t* y, jmi_int_t N) {
81    int i = 1;
82   
83    return ddot_(&N, x, &i, y, &i);
84}
85
86/* Computes  y = alpha*A*x + beta*y */
87void 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) {
88    char trans_char = trans? 'T':'N'; /* No transposition */
89    int i = 1;
90   
91    dgemv_(&trans_char, &N, &N, &a, A, &N, x, &i, &b, y, &i);
92}
93
94/* Computes A = alpha*x*y**T + A */
95void jmi_linear_algebra_dger(jmi_real_t a, jmi_real_t* x, jmi_real_t *y, jmi_real_t *A, jmi_int_t N) {
96    int i = 1;
97   
98    dger_(&N, &N, &a, x, &i, y, &i, A, &N);
99}
100
101/* Compytes z = ax + by */
102void 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) {
103    int i;
104   
105    if (a==b) {
106        for (i = 0; i < N; i++)
107            z[i] = a*(x[i]+y[i]);
108    } else if (a==-b) {
109        for (i = 0; i < N; i++)
110            z[i] = a*(x[i]-y[i]);
111    } else {
112        for (i = 0; i < N; i++)
113            z[i] = a*x[i]+b*y[i];
114    }
115}
116
117/* Find the index of the max absolute value */
118jmi_int_t jmi_linear_algebra_idamax(jmi_real_t *x, jmi_int_t N) {
119    int i = 0;
120    int j=0;
121    jmi_real_t cmax = JMI_ABS(x[i]);
122    for(i=1; i<N; i++) {
123        if(JMI_ABS(x[i])>cmax) {
124            j=i;
125        }
126    }
127    return j;
128    /*return idamax_(&N, x, &i) - 1; */ /* Compensate for Fortran indexing */
129}
130
131/* Perform LU factorization using Lapack */
132jmi_int_t jmi_linear_algebra_LU_factorize(jmi_real_t* A, jmi_int_t* pivots, jmi_int_t N) {
133    int info = 0;
134
135    dgetrf_(&N, &N, A, &N, pivots, &info);
136
137    return info;
138}
139
140/* Solve with an LU factorized matrix using Lapack */
141jmi_int_t jmi_linear_algebra_LU_solve(jmi_real_t* LU, jmi_int_t* pivots, jmi_real_t* x, jmi_int_t N) {
142    int info = 0, i = 1;
143    char trans = 'N'; /* No transposition */
144   
145    /* Back-solve and get solution in x */
146    dgetrs_(&trans, &N, &i, LU, &N, pivots, x, &N, &info);
147
148    return info;
149}
Note: See TracBrowser for help on using the repository browser.