proteus  1.8.1
C/C++/Fortran libraries
smoothers.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include "smoothers.h"
5 #include PROTEUS_LAPACK_H
6 
10 void jacobi_NR_prepare(SuperMatrix *A, double w, double tol, double* M)
11 {
12  gauss_seidel_NR_prepare(A, w, tol, M);
13 }
14 
15 void jacobi_NR_solve(SuperMatrix *A, double* M, double* R, int *node_order, double* dX)
16 {
17  /* Purpose: jacobi method with the additional assumption that A is of
18  * Stype NRformat
19  */
20  int i;
21  int N = A->nrow;
22  for (i=0; i<N; i++)
23  dX[i] = M[i]*R[i];
24 }
25 
26 void nl_jacobi_NR_solve(SuperMatrix *A, double* R, int *node_order, double w, double tol, double* dX)
27 {
28 /* Purpose: jacobi method with the additional assumption that A is of Stype NRformat
29  */
30 
31  int i;
32  int nnz;
33  double *nzval;
34  int *colind;
35  int *rowptr;
36  NRformat *AStore = (NRformat*) A->Store;
37  int N = A->nrow;
38 
39  nnz = AStore->nnz;
40  nzval = (double*)AStore->nzval;
41  colind = AStore->colind;
42  rowptr = AStore->rowptr;
43  memset(dX,0,sizeof(double)*N);
44  for (i=0; i<N; i++)
45  {
46  int cnode = node_order[i];
47  int j;
48  double diag = 1.0;
49  int diag_found = 0;
50  double val = R[cnode];
51  for (j=rowptr[cnode]; j<rowptr[cnode+1]; j++)
52  {
53  /*check for diagonal element*/
54  if (colind[j] == cnode)
55  {
56  diag = nzval[j];
57  if (fabs(nzval[j]) >= tol)
58  diag_found = 1;
59  }
60  }
61  if (diag_found)
62  dX[cnode] = w*val/diag;
63  else
64  {
65  printf("D[%d] = %12.5e \n",cnode,diag);
66  ABORT("Diagonal element is 0 within tol in Jacobi or is not in sparse matrix.");
67  }
68  }
69 }
70 
71 void gauss_seidel_NR_prepare(SuperMatrix *A, double w, double tol, double* M)
72 {
73  /* Purpose: invert diagonal for gauss_seidel method with the additional assumption that A is of Stype NRformat
74  */
75 
76  int i;
77  int nnz;
78  double *nzval;
79  int *colind;
80  int *rowptr;
81  NRformat *AStore = (NRformat*) A->Store;
82  int N = A->nrow;
83 
84  nnz = AStore->nnz;
85  nzval = (double*)AStore->nzval;
86  colind = AStore->colind;
87  rowptr = AStore->rowptr;
88 
89  for (i=0; i<N; i++)
90  {
91  int j;
92  double diag=1.0;
93  int diag_found = 0;
94  for (j=rowptr[i]; j<rowptr[i+1]; j++)
95  {
96  /*check for diagonal element*/
97  if (colind[j] == i)
98  {
99  diag = nzval[j];
100  if (fabs(nzval[j]) >= tol)
101  diag_found = 1;
102  }
103  }
104  if (diag_found)
105  M[i] = w/diag;
106  else
107  {
108  printf("D[%d] = %12.5e \n",i,diag);
109  ABORT("Diagonal element is 0 within tol in Gauss-Seidel or is not in sparse matrix.");
110  }
111  }
112 }
113 
114 void gauss_seidel_NR_solve(SuperMatrix *A, double* M, double* R, int *node_order, double* dX)
115 {
116  /* Purpose: gauss_seidel method with the additional assumption that A is of Stype NRformat
117  */
118  int i;
119  int nnz;
120  double *nzval;
121  int *colind;
122  int *rowptr;
123  NRformat *AStore = (NRformat*) A->Store;
124  int N = A->nrow;
125 
126  nnz = AStore->nnz;
127  nzval = (double*)AStore->nzval;
128  colind = AStore->colind;
129  rowptr = AStore->rowptr;
130  memset(dX,0,sizeof(double)*N);
131  for (i=0; i<N; i++)
132  {
133  int cnode = node_order[i];
134  int j;
135  double val = R[cnode];
136  for (j=rowptr[cnode]; j<rowptr[cnode+1]; j++)
137  val -= nzval[j]*dX[colind[j]];
138  dX[cnode] = M[cnode]*val;
139  }
140 }
141 
142 void nl_gauss_seidel_NR_solve(SuperMatrix *A, double* R, int *node_order, double w, double tol, double* dX)
143 {
144 /* Purpose: gauss_seidel method with the additional assumption that A is of Stype NRformat
145  */
146 
147  int i;
148  int nnz;
149  double *nzval;
150  int *colind;
151  int *rowptr;
152  NRformat *AStore = (NRformat*) A->Store;
153  int N = A->nrow;
154 
155  nnz = AStore->nnz;
156  nzval = (double*)AStore->nzval;
157  colind = AStore->colind;
158  rowptr = AStore->rowptr;
159  memset(dX,0,sizeof(double)*N);
160  for (i=0; i<N; i++)
161  {
162  int cnode = node_order[i];
163  int j;
164  double diag=1.0;
165  int diag_found = 0;
166  double val = R[cnode];
167 
168  for (j=rowptr[cnode]; j<rowptr[cnode+1]; j++)
169  {
170  /*check for diagonal element*/
171  if (colind[j] == cnode)
172  {
173  diag = nzval[j];
174  if (fabs(nzval[j]) >= tol)
175  diag_found = 1;
176  }
177  val -= nzval[j]*dX[colind[j]];
178  }
179  if (diag_found)
180  dX[cnode] = w*val/diag;
181  else
182  {
183  printf("D[%d] = %12.5e \n",cnode,diag);
184  ABORT("Diagonal element is 0 within tol in Gauss-Seidel or is not in sparse matrix.");
185  }
186  }
187 }
188 
189 int asm_NR_init(SuperMatrix *A,
190  int** subdomain_dim_p,
191  int*** l2g_L_p,
192  double*** subdomain_L_p,
193  double*** subdomain_R_p,
194  double*** subdomain_dX_p,
195  PROTEUS_LAPACK_INTEGER*** subdomain_pivots_p)
196 {
197  int i,N;
198  int* subdomain_dim;
199  int** l2g_L;
200  double** subdomain_L;
201  double** subdomain_R;
202  double** subdomain_dX;
203  PROTEUS_LAPACK_INTEGER** subdomain_pivots;
204  NRformat* ANR = (NRformat*)A->Store;
205  N = A->nrow;
206  *subdomain_dim_p = (int*)malloc(N*sizeof(int));
207  *subdomain_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
208  *l2g_L_p = (int**)malloc(N*sizeof(int*));
209  *subdomain_R_p = (double**)malloc(N*sizeof(double*));
210  *subdomain_dX_p = (double**)malloc(N*sizeof(double*));
211  *subdomain_L_p = (double**)malloc(N*sizeof(double*));
212  if ( (*subdomain_dim_p == NULL) ||
213  (*l2g_L_p == NULL) ||
214  (*subdomain_R_p == NULL) ||
215  (*subdomain_dX_p == NULL) ||
216  (*subdomain_L_p == NULL) ||
217  (*subdomain_pivots_p == NULL) )
218  {
219  return 1;
220  }
221  subdomain_dim = *subdomain_dim_p;
222  subdomain_pivots = *subdomain_pivots_p;
223  l2g_L = *l2g_L_p;
224  subdomain_R = *subdomain_R_p;
225  subdomain_dX =*subdomain_dX_p;
226  subdomain_L =*subdomain_L_p;
227  /* extract subdomain sizes and allocate storage */
228  for (i=0; i<N; i++)
229  {
230  subdomain_dim[i] = ANR->rowptr[i+1]-ANR->rowptr[i];
231  subdomain_pivots[i] = (PROTEUS_LAPACK_INTEGER*)malloc(subdomain_dim[i]*sizeof(PROTEUS_LAPACK_INTEGER));
232  l2g_L[i] = (int*)malloc(subdomain_dim[i]*subdomain_dim[i]*sizeof(int));
233  subdomain_R[i] = (double*)malloc(subdomain_dim[i]*sizeof(double));
234  subdomain_dX[i] = (double*)malloc(subdomain_dim[i]*sizeof(double));
235  subdomain_L[i] = (double*)malloc(subdomain_dim[i]*subdomain_dim[i]*sizeof(double));
236  if ((l2g_L[i] == NULL) ||
237  (subdomain_R[i] == NULL) ||
238  (subdomain_dX[i] == NULL) ||
239  (subdomain_L[i] == NULL))
240  {
241  return 1;
242  }
243  }
244  /* extract the local to global mappings */
245  for (i=0;i<N;i++)
246  {
247  int j, jj, k, kk;
248  for (kk=0; kk < subdomain_dim[i]; kk++)
249  {
250  k = ANR->colind[ANR->rowptr[i] + kk];
251  for (jj=0;jj<subdomain_dim[i]; jj++)
252  {
253  l2g_L[i][kk*subdomain_dim[i]+jj] = -1;
254  for (j=ANR->rowptr[k];j<ANR->rowptr[k+1];j++)
255  {
256  if (ANR->colind[j] == ANR->colind[ANR->rowptr[i] + jj])
257  {
258  l2g_L[i][kk*subdomain_dim[i]+jj] = j;
259  break;
260  }
261  }
262  }
263  }
264  }
265  return 0;
266 }
267 
268 void asm_NR_free(int N,
269  int* subdomain_dim,
270  int** l2g_L,
271  double** subdomain_L,
272  double** subdomain_R,
273  double** subdomain_dX,
274  PROTEUS_LAPACK_INTEGER** subdomain_pivots)
275 {
276  int i;
277  free(subdomain_dim);
278  for (i=0;i<N;i++)
279  {
280  free(subdomain_pivots[i]);
281  free(l2g_L[i]);
282  free(subdomain_R[i]);
283  free(subdomain_dX[i]);
284  free(subdomain_L[i]);
285  }
286  free(subdomain_pivots);
287  free(l2g_L);
288  free(subdomain_R);
289  free(subdomain_dX);
290  free(subdomain_L);
291 }
292 
293 void asm_NR_prepare(SuperMatrix *A,
294  int* subdomain_dim,
295  int** l2g_L,
296  double** subdomainL,
297  PROTEUS_LAPACK_INTEGER** subdomainPivots)
298 {
299  /* Purpose: extract subdomain matrices and factor in place
300  */
301  int i;
302  NRformat *ANR = (NRformat*) A->Store;
303  double *nzval = (double*) ANR->nzval;
304  for (i=0; i<A->nrow; i++)
305  {
306  int ii,jj;
307  for (ii=0;ii<subdomain_dim[i];ii++)
308  for (jj=0;jj<subdomain_dim[i];jj++)
309  {
310  int index = ii*subdomain_dim[i] + jj;
311  if (l2g_L[i][index] != -1)
312  {
313  subdomainL[i][ii + jj*subdomain_dim[i]] =
314  nzval[l2g_L[i][index]];
315  }
316  else
317  subdomainL[i][ii+jj*subdomain_dim[i]] = 0.0;
318  }
319  PROTEUS_LAPACK_INTEGER La_N=((PROTEUS_LAPACK_INTEGER)subdomain_dim[i]),INFO=0;
320  dgetrf_(&La_N,&La_N,subdomainL[i],&La_N,subdomainPivots[i],&INFO);
321  }
322 }
323 
324 void asm_NR_solve(SuperMatrix *A,
325  double w,
326  double** subdomainL,
327  int* subdomain_dim,
328  int** l2g_L,
329  double* R,
330  double** subdomainR,
331  int *node_order,
332  double** subdomain_dX,
333  double* dX,
334  PROTEUS_LAPACK_INTEGER** subdomainPivots)
335 {
336 /* Purpose: asm method with the additional assumption that A is of Stype NRformat
337  */
338  int i;
339  int nnz;
340  double *nzval;
341  int *colind;
342  int *rowptr;
343  NRformat *AStore = (NRformat*) A->Store;
344  int N = A->nrow;
345 
346  PROTEUS_LAPACK_INTEGER La_N,INFO=0,NRHS=1;
347  char TRANS='N';
348 
349  nnz = AStore->nnz;
350  nzval = (double*)AStore->nzval;
351  colind = AStore->colind;
352  rowptr = AStore->rowptr;
353  memset(dX,0,sizeof(double)*N);
354  for (i=0; i<N; i++)
355  {
356  int cnode = node_order[i];
357  int j, k,jj, ii;
358  /* extract and update the subdomain residual */
359  for (jj = 0;jj<subdomain_dim[cnode];jj++)
360  {
361  subdomainR[cnode][jj] = R[colind[rowptr[cnode] + jj]];
362  }
363  for (ii=0;ii<subdomain_dim[cnode];ii++)
364  {
365  k = colind[rowptr[cnode]+ii];
366  for (j=rowptr[k];j<rowptr[k+1];j++)
367  subdomainR[cnode][ii] -= nzval[j]*dX[colind[j]];
368  }
369  /* copy R into dX because lapack wants it that way*/
370  for (jj = 0;jj<subdomain_dim[cnode];jj++)
371  {
372  subdomain_dX[cnode][jj] = subdomainR[cnode][jj];
373  }
374  /* solve subdomain problem*/
375  La_N = (PROTEUS_LAPACK_INTEGER) subdomain_dim[cnode];
376  dgetrs_(&TRANS,
377  &La_N,
378  &NRHS,
379  subdomainL[cnode],
380  &La_N,
381  subdomainPivots[cnode],
382  subdomain_dX[cnode],
383  &La_N,
384  &INFO);
385  /* set the global correction from the subdomain correction */
386  for (jj = 0;jj<subdomain_dim[cnode];jj++)
387  dX[colind[rowptr[cnode] + jj]] += w*subdomain_dX[cnode][jj];
388  }
389 }
390 
391 int basm_NR_init(int rowBlocks,
392  SuperMatrix *A,
393  int** subdomain_dim_p,
394  int*** l2g_L_p,
395  double*** subdomain_L_p,
396  double*** subdomain_R_p,
397  double*** subdomain_dX_p,
398  PROTEUS_LAPACK_INTEGER*** subdomain_pivots_p,
399  PROTEUS_LAPACK_INTEGER*** subdomain_col_pivots_p)
400 {
401  int i,N;
402  int* subdomain_dim;
403  int** l2g_L;
404  double** subdomain_L;
405  double** subdomain_R;
406  double** subdomain_dX;
407  PROTEUS_LAPACK_INTEGER** subdomain_pivots;
408  PROTEUS_LAPACK_INTEGER** subdomain_col_pivots;
409  NRformat* ANR = (NRformat*)A->Store;
410  /*local variables for computing dimensions etc*/
411  int max_local_dim=0; /*size of largest possible subdomain system (assuming all the column entries are unique)*/
412  int local_nonzero_entries = 0; /*total number of nonzero entries in local row system*/
413  int n_unique_local_dofs=0;/*number of unique unknowns in a local
414  system*/
415  int *unique_local_dofs; /*use to keep track of which unknowns
416  have already been accounted for */
417  assert((A->nrow % rowBlocks) == 0);
418  N = A->nrow/rowBlocks;
419 
420  for (i=0; i < N; i++)
421  {
422  local_nonzero_entries = ANR->rowptr[(i+1)*rowBlocks]-ANR->rowptr[i*rowBlocks];
423  if (local_nonzero_entries > max_local_dim)
424  max_local_dim = local_nonzero_entries;
425  }
426  unique_local_dofs = (int*) malloc(max_local_dim*sizeof(int));
427  for (i=0; i < max_local_dim; i++)
428  unique_local_dofs[i]=-1;
429 
430  *subdomain_dim_p = (int*)malloc(N*sizeof(int));
431  *subdomain_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
432  *subdomain_col_pivots_p = (PROTEUS_LAPACK_INTEGER**)malloc(N*sizeof(PROTEUS_LAPACK_INTEGER*));
433  *l2g_L_p = (int**)malloc(N*sizeof(int*));
434  *subdomain_R_p = (double**)malloc(N*sizeof(double*));
435  *subdomain_dX_p = (double**)malloc(N*sizeof(double*));
436  *subdomain_L_p = (double**)malloc(N*sizeof(double*));
437  if ( (*subdomain_dim_p == NULL) ||
438  (*l2g_L_p == NULL) ||
439  (*subdomain_R_p == NULL) ||
440  (*subdomain_dX_p == NULL) ||
441  (*subdomain_L_p == NULL) ||
442  (*subdomain_pivots_p == NULL) ||
443  (*subdomain_col_pivots_p == NULL))
444  {
445  /*clean up*/
446  free(unique_local_dofs);
447  return 1;
448  }
449  subdomain_dim = *subdomain_dim_p;
450  subdomain_pivots = *subdomain_pivots_p;
451  subdomain_col_pivots = *subdomain_col_pivots_p;
452  l2g_L = *l2g_L_p;
453  subdomain_R = *subdomain_R_p;
454  subdomain_dX =*subdomain_dX_p;
455  subdomain_L =*subdomain_L_p;
456  /* extract subdomain sizes and allocate storage */
457  for (i=0; i<N; i++)
458  {
459  /*count number of unique dofs in the local system*/
460  int jj,kk,j,found_j;
461  local_nonzero_entries = ANR->rowptr[(i+1)*rowBlocks]-ANR->rowptr[i*rowBlocks];
462  for (jj=0; jj < local_nonzero_entries; jj++)
463  unique_local_dofs[jj] = -12345;
464  n_unique_local_dofs = 0;
465  for (jj=0; jj < local_nonzero_entries; jj++)
466  {
467  /*global column id*/
468  j = ANR->colind[ANR->rowptr[i*rowBlocks]+jj];
469  /*is this entry unique*/
470  found_j = -1;
471  for (kk=0; kk < n_unique_local_dofs; kk++)
472  {
473  if (unique_local_dofs[kk] == j)
474  {
475  found_j = 1;
476  break;
477  }
478  }
479  if (found_j == -1)
480  {
481  unique_local_dofs[n_unique_local_dofs] = j;
482  n_unique_local_dofs++;
483  }
484  assert(n_unique_local_dofs <= local_nonzero_entries);
485  }
486 
487  subdomain_dim[i] = n_unique_local_dofs;
488  subdomain_pivots[i] = (PROTEUS_LAPACK_INTEGER*)malloc(subdomain_dim[i]*sizeof(PROTEUS_LAPACK_INTEGER));
489  subdomain_col_pivots[i] = (PROTEUS_LAPACK_INTEGER*)malloc(subdomain_dim[i]*sizeof(PROTEUS_LAPACK_INTEGER));
490  l2g_L[i] = (int*)malloc(subdomain_dim[i]*subdomain_dim[i]*sizeof(int));
491  subdomain_R[i] = (double*)malloc(subdomain_dim[i]*sizeof(double));
492  subdomain_dX[i] = (double*)malloc(subdomain_dim[i]*sizeof(double));
493  subdomain_L[i] = (double*)malloc(subdomain_dim[i]*subdomain_dim[i]*sizeof(double));
494  if ((l2g_L[i] == NULL) ||
495  (subdomain_R[i] == NULL) ||
496  (subdomain_dX[i] == NULL) ||
497  (subdomain_L[i] == NULL))
498  {
499  /*clean up*/
500  free(unique_local_dofs);
501  return 1;
502  }
503  }
504  /* extract the local to global mappings */
505  for (i=0;i<N;i++)
506  {
507  int j, jj, k, kk, found_j;
508  /*again determine unique column indeces*/
509  local_nonzero_entries = ANR->rowptr[(i+1)*rowBlocks]-ANR->rowptr[i*rowBlocks];
510  for (jj=0; jj < local_nonzero_entries; jj++)
511  unique_local_dofs[jj] = -12345;
512  n_unique_local_dofs = 0;
513  for (jj=0; jj < local_nonzero_entries; jj++)
514  {
515  /*global column id*/
516  j = ANR->colind[ANR->rowptr[i*rowBlocks]+jj];
517  /*is this entry unique*/
518  found_j = -1;
519  for (kk=0; kk < n_unique_local_dofs; kk++)
520  {
521  if (unique_local_dofs[kk] == j)
522  {
523  found_j = 1;
524  break;
525  }
526  }
527  if (found_j == -1)
528  {
529  unique_local_dofs[n_unique_local_dofs] = j;
530  n_unique_local_dofs++;
531  }
532  assert(n_unique_local_dofs <= local_nonzero_entries);
533  }
534  assert(subdomain_dim[i] == n_unique_local_dofs);
535  for (kk=0; kk < subdomain_dim[i]; kk++)
536  {
537  k = unique_local_dofs[kk];
538  for (jj=0;jj<subdomain_dim[i]; jj++)
539  {
540  l2g_L[i][kk*subdomain_dim[i]+jj] = -1;
541  for (j=ANR->rowptr[k];j<ANR->rowptr[k+1];j++)/*just look through k's row*/
542  {
543  if (ANR->colind[j] == unique_local_dofs[jj])
544  {
545  l2g_L[i][kk*subdomain_dim[i]+jj] = j;
546  break;
547  }
548  }
549  }
550  }
551  }
552  /*mwf debug*/
553  printf("BASM_init N= %d rowBlocks=%d \n",N,rowBlocks);
554  for (i=0; i < N; i++)
555  {
556  int kk,jj;
557  printf("l2g_L[%d] dim= %d\n",i,subdomain_dim[i]);
558  for (kk=0; kk < subdomain_dim[i]; kk++)
559  {
560  printf("\t");
561  for (jj=0; jj < subdomain_dim[i]; jj++)
562  {
563  printf(" %d ",l2g_L[i][kk*subdomain_dim[i] + jj]);
564  }
565  printf("\n");
566  }
567  }
568  /*clean up*/
569  free(unique_local_dofs);
570 
571  return 0;
572 }
573 
574 void basm_NR_free(int N,
575  int* subdomain_dim,
576  int** l2g_L,
577  double** subdomain_L,
578  double** subdomain_R,
579  double** subdomain_dX,
580  PROTEUS_LAPACK_INTEGER** subdomain_pivots,
581  PROTEUS_LAPACK_INTEGER** subdomain_col_pivots)
582 {
583  int i;
584  free(subdomain_dim);
585  for (i=0;i<N;i++)
586  {
587  free(subdomain_pivots[i]);
588  free(subdomain_col_pivots[i]);
589  free(l2g_L[i]);
590  free(subdomain_R[i]);
591  free(subdomain_dX[i]);
592  free(subdomain_L[i]);
593  }
594  free(subdomain_pivots);
595  free(subdomain_col_pivots);
596  free(l2g_L);
597  free(subdomain_R);
598  free(subdomain_dX);
599  free(subdomain_L);
600 }
601 
602 void basm_NR_prepare(int rowBlocks,
603  int N,
604  SuperMatrix *A,
605  int* subdomain_dim,
606  int** l2g_L,
607  double** subdomainL,
608  PROTEUS_LAPACK_INTEGER** subdomainPivots,
609  PROTEUS_LAPACK_INTEGER** subdomainColPivots)
610 {
611  /* Purpose: extract subdomain matrices and factor in place
612  */
613  int i;
614  NRformat *ANR = (NRformat*) A->Store;
615  double *nzval = (double*) ANR->nzval;
616  assert (N*rowBlocks == A->nrow);
617  /*zero things for safety*/
618  for (i=0; i< N; i++)
619  {
620  int ii,jj;
621  for (ii=0; ii<subdomain_dim[i]; ii++)
622  {
623  subdomainPivots[i][ii] = 0; subdomainColPivots[i][ii] = 0;
624  for (jj=0; jj<subdomain_dim[i]; jj++)
625  subdomainL[i][jj*subdomain_dim[i]+ii] = 0.0;
626  }
627  }
628  for (i=0; i<N; i++)
629  {
630  int ii,jj;
631  for (ii=0;ii<subdomain_dim[i];ii++)
632  for (jj=0;jj<subdomain_dim[i];jj++)
633  {
634  int index = ii*subdomain_dim[i] + jj;
635  if (l2g_L[i][index] != -1)
636  {
637  subdomainL[i][ii + jj*subdomain_dim[i]] =
638  nzval[l2g_L[i][index]];
639  }
640  else
641  subdomainL[i][ii+jj*subdomain_dim[i]] = 0.0;
642  }
643  PROTEUS_LAPACK_INTEGER La_N=((PROTEUS_LAPACK_INTEGER)subdomain_dim[i]),INFO=0;
644  dgetc2_(&La_N,subdomainL[i],&La_N,subdomainPivots[i],subdomainColPivots[i],&INFO);
645  /*doesn't seem to be handling trivial case of dim=1 and L_00 = 0 well*/
646  /*mwf debug*/
647  if (INFO > 0)
648  {
649  printf("basm prepare jac dgetc2 INFO=%d nN=%d \n",(int)(INFO),subdomain_dim[i]);
650  for (ii=0;ii<subdomain_dim[i];ii++)
651  {
652  for(jj=0;jj<subdomain_dim[i];jj++)
653  {
654 
655  printf("%12.5e \t",subdomainL[i][ii*subdomain_dim[i] + jj]);
656  }
657  printf("\n");
658  }
659  for (ii=0; ii<subdomain_dim[i];ii++)
660  {
661  printf("colPivot[%d]= %dl \t",ii,subdomainColPivots[i][ii]);
662  printf("pivot[%d]= %dl \t",ii,subdomainPivots[i][ii]);
663  }
664  }
665  }
666 }
667 
668 void basm_NR_solve(int rowBlocks,
669  int N,
670  SuperMatrix *A,
671  double w,
672  double** subdomainL,
673  int* subdomain_dim,
674  int** l2g_L,
675  double* R,
676  double** subdomainR,
677  int *node_order,
678  double** subdomain_dX,
679  double* dX,
680  PROTEUS_LAPACK_INTEGER** subdomainPivots,
681  PROTEUS_LAPACK_INTEGER** subdomainColPivots)
682 {
683 /* Purpose: asm method with the additional assumption that A is of Stype NRformat
684  */
685  int i;
686  int nnz;
687  double *nzval;
688  int *colind;
689  int *rowptr;
690  NRformat *AStore = (NRformat*) A->Store;
691  assert(N*rowBlocks == A->nrow);
692  double scale = 1.0;
693  PROTEUS_LAPACK_INTEGER La_N;
694  /*have to keep track of unique dofs in local system now*/
695  int max_local_dim=0; /*size of largest possible subdomain system (assuming all the column entries are unique)*/
696  int local_nonzero_entries = 0; /*total number of nonzero entries in local row system*/
697  int n_unique_local_dofs=0;/*number of unique unknowns in a local
698  system*/
699  int *unique_local_dofs; /*use to keep track of which unknowns
700  have already been accounted for */
701 
702 
703  for (i=0; i < N; i++)
704  {
705  local_nonzero_entries = AStore->rowptr[(i+1)*rowBlocks]-AStore->rowptr[i*rowBlocks];
706  if (local_nonzero_entries > max_local_dim)
707  max_local_dim = local_nonzero_entries;
708  }
709  unique_local_dofs = (int*) malloc(max_local_dim*sizeof(int));
710  for (i=0; i < max_local_dim; i++)
711  unique_local_dofs[i]=-1;
712 
713  nnz = AStore->nnz;
714  nzval = (double*)AStore->nzval;
715  colind = AStore->colind;
716  rowptr = AStore->rowptr;
717  memset(dX,0,sizeof(double)*N*rowBlocks);
718  for (i=0; i<N; i++)
719  {
720  int cnode = node_order[i];
721  int j, k,kk, jj, ii, found_j=-1;
722  /*again determine unique column indeces maybe go ahead and store these*/
723  local_nonzero_entries = AStore->rowptr[(cnode+1)*rowBlocks]-AStore->rowptr[cnode*rowBlocks];
724  for (jj=0; jj < local_nonzero_entries; jj++)
725  unique_local_dofs[jj] = -12345;
726  n_unique_local_dofs = 0;
727  for (jj=0; jj < local_nonzero_entries; jj++)
728  {
729  /*global column id*/
730  j = AStore->colind[AStore->rowptr[cnode*rowBlocks]+jj];
731  /*is this entry unique*/
732  found_j = -1;
733  for (kk=0; kk < n_unique_local_dofs; kk++)
734  {
735  if (unique_local_dofs[kk] == j)
736  {
737  found_j = 1;
738  break;
739  }
740  }
741  if (found_j == -1)
742  {
743  unique_local_dofs[n_unique_local_dofs] = j;
744  n_unique_local_dofs++;
745  }
746  assert(n_unique_local_dofs <= local_nonzero_entries);
747  }
748  assert(subdomain_dim[cnode] == n_unique_local_dofs);
749  /* extract and update the subdomain residual */
750  for (jj = 0;jj<subdomain_dim[cnode];jj++)
751  {
752  subdomainR[cnode][jj] = R[unique_local_dofs[jj]]; /*colind[rowptr[cnode*rowBlocks] + jj]];*/
753  }
754  for (ii=0;ii<subdomain_dim[cnode];ii++)
755  {
756  k = unique_local_dofs[ii]; /*colind[rowptr[cnode*rowBlocks]+ii];*/
757  for (j=rowptr[k];j<rowptr[k+1];j++)
758  subdomainR[cnode][ii] -= nzval[j]*dX[colind[j]];/*check this if needs to be only entries in nodestar*/
759  }
760  /* copy R into dX because lapack wants it that way*/
761  for (jj = 0;jj<subdomain_dim[cnode];jj++)
762  {
763  subdomain_dX[cnode][jj] = subdomainR[cnode][jj];
764  }
765  /* solve subdomain problem*/
766  La_N = (PROTEUS_LAPACK_INTEGER) subdomain_dim[cnode];
767  /*doesn't seem to be handling trivial case of dim=1 and L_00 = 0 well*/
768  if (La_N == 1 && fabs(subdomainL[cnode][0]) <= 1.0e-64)
769  subdomain_dX[cnode][0] = 0.0;
770  else
771  {
772  dgesc2_(&La_N,
773  subdomainL[cnode],
774  &La_N,
775  subdomain_dX[cnode],
776  subdomainPivots[cnode],
777  subdomainColPivots[cnode],
778  &scale);
779  }
780  /* set the global correction from the subdomain correction */
781  for (jj = 0;jj<subdomain_dim[cnode];jj++)
782  {
783  /*dX[colind[rowptr[cnode*rowBlocks] + jj]] += w*subdomain_dX[cnode][jj];*/
784  dX[unique_local_dofs[jj]] += w*subdomain_dX[cnode][jj];
785  }
786  }
787 }
788 
asm_NR_prepare
void asm_NR_prepare(SuperMatrix *A, int *subdomain_dim, int **l2g_L, double **subdomainL, PROTEUS_LAPACK_INTEGER **subdomainPivots)
Definition: smoothers.c:293
jacobi_NR_prepare
void jacobi_NR_prepare(SuperMatrix *A, double w, double tol, double *M)
Definition: smoothers.c:10
w
#define w(x)
Definition: jf.h:22
basm_NR_prepare
void basm_NR_prepare(int rowBlocks, int N, SuperMatrix *A, int *subdomain_dim, int **l2g_L, double **subdomainL, PROTEUS_LAPACK_INTEGER **subdomainPivots, PROTEUS_LAPACK_INTEGER **subdomainColPivots)
Definition: smoothers.c:602
smoothers.h
C implementations of multilevel smoother algorithms.
jacobi_NR_solve
void jacobi_NR_solve(SuperMatrix *A, double *M, double *R, int *node_order, double *dX)
Definition: smoothers.c:15
gauss_seidel_NR_prepare
void gauss_seidel_NR_prepare(SuperMatrix *A, double w, double tol, double *M)
Definition: smoothers.c:71
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
asm_NR_free
void asm_NR_free(int N, int *subdomain_dim, int **l2g_L, double **subdomain_L, double **subdomain_R, double **subdomain_dX, PROTEUS_LAPACK_INTEGER **subdomain_pivots)
Definition: smoothers.c:268
dgetrf_
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info)
R
Double R
Definition: Headers.h:82
basm_NR_init
int basm_NR_init(int rowBlocks, SuperMatrix *A, int **subdomain_dim_p, int ***l2g_L_p, double ***subdomain_L_p, double ***subdomain_R_p, double ***subdomain_dX_p, PROTEUS_LAPACK_INTEGER ***subdomain_pivots_p, PROTEUS_LAPACK_INTEGER ***subdomain_col_pivots_p)
Definition: smoothers.c:391
asm_NR_init
int asm_NR_init(SuperMatrix *A, int **subdomain_dim_p, int ***l2g_L_p, double ***subdomain_L_p, double ***subdomain_R_p, double ***subdomain_dX_p, PROTEUS_LAPACK_INTEGER ***subdomain_pivots_p)
Definition: smoothers.c:189
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
asm_NR_solve
void asm_NR_solve(SuperMatrix *A, double w, double **subdomainL, int *subdomain_dim, int **l2g_L, double *R, double **subdomainR, int *node_order, double **subdomain_dX, double *dX, PROTEUS_LAPACK_INTEGER **subdomainPivots)
Definition: smoothers.c:324
basm_NR_free
void basm_NR_free(int N, int *subdomain_dim, int **l2g_L, double **subdomain_L, double **subdomain_R, double **subdomain_dX, PROTEUS_LAPACK_INTEGER **subdomain_pivots, PROTEUS_LAPACK_INTEGER **subdomain_col_pivots)
Definition: smoothers.c:574
nl_gauss_seidel_NR_solve
void nl_gauss_seidel_NR_solve(SuperMatrix *A, double *R, int *node_order, double w, double tol, double *dX)
Definition: smoothers.c:142
dgetrs_
int dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info)
nl_jacobi_NR_solve
void nl_jacobi_NR_solve(SuperMatrix *A, double *R, int *node_order, double w, double tol, double *dX)
Definition: smoothers.c:26
gauss_seidel_NR_solve
void gauss_seidel_NR_solve(SuperMatrix *A, double *M, double *R, int *node_order, double *dX)
Definition: smoothers.c:114
basm_NR_solve
void basm_NR_solve(int rowBlocks, int N, SuperMatrix *A, double w, double **subdomainL, int *subdomain_dim, int **l2g_L, double *R, double **subdomainR, int *node_order, double **subdomain_dX, double *dX, PROTEUS_LAPACK_INTEGER **subdomainPivots, PROTEUS_LAPACK_INTEGER **subdomainColPivots)
Definition: smoothers.c:668
nnz
#define nnz
Definition: Richards.h:9
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20