1/*
2	luserver.h
3
4	LU factorization Web service.
5
6--------------------------------------------------------------------------------
7gSOAP XML Web services tools
8Copyright (C) 2001-2008, Robert van Engelen, Genivia, Inc. All Rights Reserved.
9This software is released under one of the following two licenses:
10GPL or Genivia's license for commercial use.
11--------------------------------------------------------------------------------
12GPL license.
13
14This program is free software; you can redistribute it and/or modify it under
15the terms of the GNU General Public License as published by the Free Software
16Foundation; either version 2 of the License, or (at your option) any later
17version.
18
19This program is distributed in the hope that it will be useful, but WITHOUT ANY
20WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
21PARTICULAR PURPOSE. See the GNU General Public License for more details.
22
23You should have received a copy of the GNU General Public License along with
24this program; if not, write to the Free Software Foundation, Inc., 59 Temple
25Place, Suite 330, Boston, MA 02111-1307 USA
26
27Author contact information:
28engelen@genivia.com / engelen@acm.org
29--------------------------------------------------------------------------------
30A commercial use license is available from Genivia, Inc., contact@genivia.com
31--------------------------------------------------------------------------------
32*/
33
34#include "soapH.h"
35#include <math.h>
36
37////////////////////////////////////////////////////////////////////////////////
38//
39//  main
40//
41////////////////////////////////////////////////////////////////////////////////
42
43int main(int argc, char **argv)
44{ struct soap *soap;
45  int m, s;
46  soap = soap_new();
47  if (argc < 3)
48    soap_serve(soap); // run as CGI application over the Web
49  else // run as stand-alone server on machine given by argv[1] listening to port argv[2]
50  { m = soap_bind(soap, argv[1], atoi(argv[2]), 100);
51    if (m < 0)
52    { soap_print_fault(soap, stderr);
53      exit(-1);
54    }
55    fprintf(stderr, "Socket connection successful: master socket = %d\n", m);
56    for (int i = 1; ; i++)
57    { s = soap_accept(soap);
58      if (s < 0)
59      { soap_print_fault(soap, stderr);
60        exit(-1);
61      }
62      fprintf(stderr, "%d: accepted connection from IP = %d.%d.%d.%d socket = %d ... ", i, (int)(soap->ip>>24)&0xFF, (int)(soap->ip>>16)&0xFF, (int)(soap->ip>>8)&0xFF, (int)soap->ip&0xFF, s);
63      soap_serve(soap);		// process request
64      fprintf(stderr, "request served\n");
65      soap_destroy(soap);	// delete class instances
66      soap_end(soap);		// clean up everything and close socket
67    }
68  }
69  soap_done(soap);
70  return 0;
71}
72
73////////////////////////////////////////////////////////////////////////////////
74//
75//  LU decomposition: remote method interface
76//
77////////////////////////////////////////////////////////////////////////////////
78
79int ludcmp(struct soap*, matrix&, ivector&, double&);
80
81int ns1__ludcmp(struct soap *soap, matrix *a, struct ns1__ludcmpResponse &result)
82{ result.a = a;
83  result.i = soap_new_ivector(soap, -1);
84  if (ludcmp(soap, *result.a, *result.i, result.d))
85    return soap_sender_fault(soap, "Singular matrix in routine ludcmp", NULL);
86  return SOAP_OK;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90//
91//  LU decomposition: algorithm
92//
93////////////////////////////////////////////////////////////////////////////////
94
95int ludcmp(struct soap *soap, matrix &a, ivector &indx, double &d)
96{ int i, imax = 0, j, k, n;
97  double big, dum, sum, temp;
98  n = a.size();
99  vector vv(soap);
100  vv.resize(n);
101  indx.resize(n);
102  d = 1.0;
103  for (i = 1; i <= n; i++)
104  { big = 0.0;
105    a[i].resize(n);
106    for (j = 1; j <= n; j++)
107      if ((temp = fabs(a[i][j])) > big)
108        big = temp;
109    if (big == 0.0)
110      return -1;
111    vv[i] = 1.0/big;
112  }
113  for (j = 1; j <= n; j++)
114  { for (i = 1; i < j; i++)
115    { sum = a[i][j];
116      for (k = 1; k < i; k++)
117        sum -= a[i][k]*a[k][j];
118      a[i][j] = sum;
119    }
120    big = 0.0;
121    for (i = j; i <= n; i++)
122    { sum = a[i][j];
123      for (k = 1; k < j; k++)
124        sum -= a[i][k]*a[k][j];
125      a[i][j] = sum;
126      if ((dum = vv[i]*fabs(sum)) >= big)
127      { big = dum;
128        imax = i;
129      }
130    }
131    if (j != imax)
132    { for (k = 1; k <= n; k++)
133      { dum = a[imax][k];
134        a[imax][k] = a[j][k];
135        a[j][k] = dum;
136      }
137      d = -d;
138      vv[imax] = vv[j];
139    }
140    indx[j] = imax;
141    if (a[j][j] == 0.0)
142      a[j][j] = 1.0e-20;
143    if (j != n)
144    { dum = 1.0/a[j][j];
145      for (i = j+1; i <= n; i++)
146        a[i][j] *= dum;
147    }
148  }
149  for (i = 1; i <= n; i++)
150  { for (j = 1; j <= n; j++)
151      if (fabs(a[i][j]) > 1.0e-15)
152        break;
153    for (k = n; k > j; k--)
154      if (fabs(a[i][k]) > 1.0e-15)
155        break;
156    a[i].resize(j, k);
157  }
158  return 0;
159}
160
161////////////////////////////////////////////////////////////////////////////////
162//
163//  Forward- and backsubstitution: remote method interface
164//
165////////////////////////////////////////////////////////////////////////////////
166
167int lubksb(matrix&, ivector&, vector &b);
168
169int ns1__lubksb(struct soap *soap, matrix *a, ivector *i, vector *b, vector *x)
170{ lubksb(*a, *i, *b);
171  *x = *b;
172  return SOAP_OK;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176//
177//  Forward- and backsubstitution: algorithm
178//
179////////////////////////////////////////////////////////////////////////////////
180
181int lubksb(matrix &a, ivector &indx, vector &b)
182{ int i, j, k, ip, n, m, ii = 0;
183  double sum;
184  n = a.size();
185  b.resize(n);
186  for (i = 1; i <= n; i++)
187  { ip = indx[i];
188    sum = b[ip];
189    b[ip] = b[i];
190    if (ii)
191    { k = a[i].start();
192      if (ii > k)
193        k = ii;
194      m = a[i].end();
195      if (i-1 < m)
196        m = i-1;
197      for (j = k; j <= m; j++)
198        sum -= a[i][j]*b[j];
199    }
200    else if (sum)
201      ii = i;
202    b[i] = sum;
203  }
204  for (i = n; i >= 1; i--)
205  { sum = b[i];
206    k = a[i].start();
207    if (i+1 > k)
208      k = i+1;
209    m = a[i].end();
210    if (n < m)
211      m = n;
212    for (j = k; j <= m; j++)
213      sum -= a[i][j]*b[j];
214    b[i] = sum/a[i][i];
215  }
216  return SOAP_OK;
217}
218
219////////////////////////////////////////////////////////////////////////////////
220//
221//  LU solve: remote method interface
222//
223////////////////////////////////////////////////////////////////////////////////
224
225int ns1__lusol(struct soap *soap, matrix *a, vector *b, vector *x)
226{ ivector indx(soap);
227  double d;
228  if (ludcmp(soap, *a, indx, d))
229    return soap_sender_fault(soap, "Singular matrix in routine ludcmp", NULL);
230  lubksb(*a, indx, *b);
231  *x = *b;
232  return SOAP_OK;
233}
234
235////////////////////////////////////////////////////////////////////////////////
236//
237//  LU solve multiple: remote method interface
238//
239////////////////////////////////////////////////////////////////////////////////
240
241int ns1__lusols(struct soap *soap, matrix *a, matrix *b, matrix *x)
242{ ivector indx(soap);
243  double d;
244  if (ludcmp(soap, *a, indx, d))
245    return soap_sender_fault(soap, "Singular matrix in routine ludcmp", NULL);
246  for (int i = 1; i <= b->size(); i++)
247    lubksb(*a, indx, (*b)[i]);
248  *x = *b;
249  return SOAP_OK;
250}
251
252////////////////////////////////////////////////////////////////////////////////
253//
254//  LU inverse: remote method interface
255//
256////////////////////////////////////////////////////////////////////////////////
257
258int ns1__luinv(struct soap *soap, matrix *a, matrix *b)
259{ vector col(soap);
260  ivector indx(soap);
261  double d;
262  int i, j, k, n;
263  if (ludcmp(soap, *a, indx, d))
264    return soap_sender_fault(soap, "Singular matrix in routine ludcmp", NULL);
265  n = a->size();
266  col.resize(n);
267  b->resize(n, n);
268  for (j = 1; j <= n; j++)
269  { for (i = 1; i <= n; i++)
270      col[i] = 0.0;
271    col[j] = 1.0;
272    lubksb(*a, indx, col);
273    for (i = 1; i <= n; i++)
274      (*b)[i][j] = col[i];
275  }
276  for (i = 1; i <= n; i++)
277  { for (j = 1; j <= n; j++)
278      if (fabs((*b)[i][j]) > 1.0e-15)
279        break;
280    for (k = n; k > j; k--)
281      if (fabs((*b)[i][k]) > 1.0e-15)
282        break;
283    (*b)[i].resize(j, k);
284  }
285  return SOAP_OK;
286}
287
288////////////////////////////////////////////////////////////////////////////////
289//
290//  LU determinant: remote method interface
291//
292////////////////////////////////////////////////////////////////////////////////
293
294int ns1__ludet(struct soap *soap, matrix *a, double &d)
295{ ivector indx(soap);
296  if (ludcmp(soap, *a, indx, d))
297    return soap_sender_fault(soap, "Singular matrix in routine ludcmp", NULL);
298  for (int i = 1; i <= a->__size; i++)
299    d *= (*a)[i][i];
300  return SOAP_OK;
301}
302
303struct Namespace namespaces[] =
304{// "ns-prefix", "ns-name", "ns-pattern"
305  {"SOAP-ENV", "http://schemas.xmlsoap.org/soap/envelope/"},
306  {"SOAP-ENC", "http://schemas.xmlsoap.org/soap/encoding/"},
307  {"xsi", "http://www.w3.org/2001/XMLSchema-instance", "http://wwww.w3.org/*/XMLSchema-instance"},
308  {"xsd", "http://www.w3.org/2001/XMLSchema", "http://www.w3.org/*/XMLSchema"},
309  {"ns1", "urn:lu"},
310  {NULL, NULL}
311};
312