proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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
10void jacobi_NR_prepare(SuperMatrix *A, double w, double tol, double* M)
11{
12 gauss_seidel_NR_prepare(A, w, tol, M);
13}
14
15void 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
26void 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
71void 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
114void 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
142void 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
189int 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
268void 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
293void 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
324void 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
391int 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
574void 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
602void 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
668void 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
R
Double R
Definition Headers.h:82
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
gauss_seidel_NR_prepare
void gauss_seidel_NR_prepare(SuperMatrix *A, double w, double tol, double *M)
Definition smoothers.c:71
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
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
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_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
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_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
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
jacobi_NR_prepare
void jacobi_NR_prepare(SuperMatrix *A, double w, double tol, double *M)
Definition smoothers.c:10
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
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
jacobi_NR_solve
void jacobi_NR_solve(SuperMatrix *A, double *M, double *R, int *node_order, double *dX)
Definition smoothers.c:15
w
#define w(x)
Definition jf.h:22
nnz
#define nnz
Definition m_comp_co2.h:19
dgetrf_
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info)
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
dgetrs_
int dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info)
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
smoothers.h
C implementations of multilevel smoother algorithms.
Generated on for proteus by doxygen 1.18.0