TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
linear_span.c
Go to the documentation of this file.
1 
2 /* MATH2 CONTAINS ALGEBRAIC ROUTINES FOR GEOtop AND OTHER MODELS
3 MATH2 Version 0.9375 KMackenzie
4 
5 file linear_span.c
6 
7 Copyright, 2009 Stefano Endrizzi, Emanuele Cordano, Matteo Dall'Amico and Riccardo Rigon
8 
9 This file is part of MATH2.
10  MATH2 is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  MATH2 is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 
25 
26 
27 #include "turtle.h"
28 #include "tensor3D.h"
29 //#include "networks.h"
30 #include "t_utilities.h"
31 #include "linear_span.h"
32 //#include "gridded.element.input.geotop.h"
33 #include "rw_maps.h"
34 
35 
36 int no_value_function(double x,DOUBLEVECTOR *V) {
51  if ((V->element[1]>0 && x>=V->element[2]) || (V->element[1]<0 && x<=V->element[2])) {
53  return 1;
54  }else{
56  return 0;
57  }
58 
59 }
60 
61 DOUBLEMATRIX *extract_a_new_map(DOUBLETENSOR *xtensor, long l,T_INIT *UVref) {
80  long r,c;
81  DOUBLEMATRIX *xm;
82  char *error_message;
83  xm=new_doublematrix(xtensor->nrh,xtensor->nch);
84  if (l<xtensor->ndl && l>xtensor->ndh) {
85 
86  sprintf(error_message,"Number of layers does not correspond in extract_new_map function with doubletensor NAME: %s l=%ld ndl=%ld ndh=%ld \n",xtensor->name,l,xtensor->ndl,xtensor->ndh);
87  t_error(error_message);
88 
89  }else{
90 
91  for(r=xm->nrl;r<=xm->nrh;r++){
92  for(c=xm->ncl;c<=xm->nch;c++){
93  if (no_value_function(xtensor->element[l][r][c],UVref->V)==1){
94  xm->element[r][c]=UVref->V->element[2];
95  }else{
96  xm->element[r][c]=xtensor->element[l][r][c];
97  }
98  }
99 
100  }
101  return xm;
102  }
103 
104  return NULL;
105  }
106 
107 
108 
109 DOUBLEVECTOR *prod_doublematvet(DOUBLEMATRIX *m, DOUBLEVECTOR *v)
110 {
111 
112 
122  DOUBLEVECTOR *p;
123  p=new_doublevector (m->nch);
124  int i,j;
125  if (m->nch != v->nh) t_error("Error in prod_doublematvet(): The matrix and the vector have not proper dimensions\n");
126  for(i=m->nrl;i<=m->nrh;i++) {
127  double buf=0.0;
128  for (j=m->ncl;j<=m->nch;j++) {
129  buf+=m->element[i][j]*v->element[j];
130  p->element[i]=buf;
131  }
132  }
133  return p;
134 }
135 
136 
137 double prodscal (DOUBLEVECTOR *a, DOUBLEVECTOR *b)
150 {
151  double p=0.0;
152  long i;
153  if (a->nh != b->nh)
154  t_error("Error in prodscal(): The two vectors have not equal dimensions\n");
155  for(i=a->nl;i<=a->nh;i++) {
156  p+=a->element[i]*b->element[i];
157  }
158  return p;
159 }
160 
161 DOUBLEVECTOR *scalxvet (double a, DOUBLEVECTOR *b)
163 {
164  DOUBLEVECTOR *p;
165  p=new_doublevector(b->nh);
166  int i;
167  for(i=b->nl;i<=b->nh;i++) {
168  p->element[i]= a * b->element[i];
169  }
170  return p;
171 }
172 
194 DOUBLETENSOR *linear_span_doubletensor(double c1, double c2, DOUBLETENSOR *T1,DOUBLETENSOR *T2, DOUBLEVECTOR *V){
195 
196  long r,c,l;
197  DOUBLETENSOR* TL;
198 
199  if (T1->nrh!=T2->nrh) {
200  t_error("Error in function linear_span_doubletensor: doubltensors do not have the same numbers of rows!! ");
201  }else if (T1->nch!=T2->nch){
202  t_error("Error in function linear_span_doubletensor: doubletensors do not have the same numbers of columns!! ");
203  }else if(T1->ndh!=T2->ndh) {
204  t_error("linear_span_doubletensor: doubletensors do not have the same numbers of columns!!");
205  } else{
206  TL=new_doubletensor(T1->ndh,T1->nrh,T1->nch);
207  for(r=T1->nrl;r<=T1->nrh;r++){
208  for(c=T1->ncl;c<=T1->nch;c++){
209  for(l=T1->ndl;l<=T1->ndh;l++){
210  if (no_value_function(T1->element[l][r][c],V)==1 && no_value_function(T2->element[l][r][c],V)==1){
211  TL->element[l][r][c]=V->element[2];
212  }else {
213  TL->element[l][r][c]=V->element[2]=c1*T1->element[l][r][c]+c2*T2->element[l][r][c];
214  }
215  }
216  }
217 
218  }
219  return TL;
220  }
221 return NULL;
222 
223 }
224 
225 
226 DOUBLEMATRIX *linear_span_doublematrix(double c1, double c2, DOUBLEMATRIX *M1a,DOUBLEMATRIX *M2a, DOUBLEVECTOR *V){
249  long r,c;
250  DOUBLEMATRIX *ML;
251 
252  if (M1a->nrh!=M2a->nrh) {
253  t_error("Error in function linear_span_doublematrix: doublematrices do not have the same numbers of rows!! ");
254  }else if (M1a->nch!=M2a->nch) t_error("Error in function linear_span_doublematrix: doublematrices do not have the same numbers of columns!! ");
255 
256  ML=new_doublematrix(M1a->nrh,M1a->nch);
257  for(r=M1a->nrl;r<=M1a->nrh;r++){
258  for(c=M1a->ncl;c<=M1a->nch;c++){
259  if (no_value_function(M1a->element[r][c],V)==1 && no_value_function(M2a->element[r][c],V)==1){
260  ML->element[r][c]=V->element[2];
261  }else {
262  ML->element[r][c]=V->element[2]=c1*M1a->element[r][c]+c2*M2a->element[r][c];
263 
264  }
265  }
266 
267  }
268  return ML;
269  }
270 
271 
272 DOUBLEMATRIX *transpose_doublematrix(DOUBLEMATRIX *M){
285  DOUBLEMATRIX *MR;
286  long r,c;
287 
288  MR=new_doublematrix(M->nch,M->nrh);
289  for(r=MR->nrl;r<=MR->nrh;r++){
290  for(r=MR->nrl;r<=MR->nrh;r++){
291  MR->element[r][c]=M->element[c][r];
292  }
293 
294  }
295 
296  return MR;
297 
298 }
299 
300 DOUBLEVECTOR *extract_a_column_from_doublematrix(long d,DOUBLEMATRIX *M) {
315  DOUBLEVECTOR *VC;
316  long r;
317 
318  VC=new_doublevector(M->nrh);
319  if ((d>=M->ncl) && (d<=M->nch)) {
320  for(r=M->nrl;r<=M->nrh;r++){
321  VC->element[r]=M->element[r][d];
322  }
323  } else {
324  printf ("Warning: d=%ld is out of range: %ld to %ld columns",d,M->ncl,M->nch);
325  }
326 
327  return VC;
328 
329 
330 }
331 
332 
333 DOUBLEVECTOR *extract_a_row_from_doublematrix(long d,DOUBLEMATRIX *M) {
348  DOUBLEVECTOR *VC;
349  long c;
350 
351  VC=new_doublevector(M->nch);
352  if ((d>=M->nrl) && (d<=M->nrh)) {
353  for(c=M->nrl;c<=M->nrh;c++){
354  VC->element[c]=M->element[d][c];
355  }
356  } else {
357  printf ("Warning: d=%ld is out of range: %ld to %ld rows",d,M->nrl,M->nrh);
358  }
359 
360  return VC;
361 
362 
363 }
364 
365 
366 
367 
368 DOUBLEVECTOR *extract_a_vertical_column_from_doubletensor(long r,long c, DOUBLETENSOR *T) {
383  DOUBLEVECTOR *VC;
384  long i;
385 
386  VC=new_doublevector(T->ndh);
387  if ((c>=T->ncl) && (c<=T->nch) && (r>=T->nrl) && (r<=T->nrh) ) {
388  for(i=VC->nl;i<=VC->nh;i++){
389  VC->element[i]=T->element[i][r][c];
390  }
391  } else {
392  printf ("Warning: r=%ld c=%ld is out of range: %ld to %ld columns and %ld to %ld rows ",r,c,T->nrl,T->nrh,T->ncl,T->nch);
393  }
394 
395  return VC;
396 
397 
398 
399 }
400 
401 
402 
403 
404