home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / sparseio < prev    next >
Encoding:
Text File  |  1994-01-13  |  8.0 KB  |  340 lines
  1. >elt[idx1].nxt_row = tmp_e->nxt_row;
  2.         r1->elt[idx1].nxt_idx = tmp_e->nxt_idx;
  3.         tmp_e->nxt_row = i1;
  4.         tmp_e->nxt_idx = idx1;
  5.     }
  6.     }
  7.     else if ( r1->elt[idx1].col != j1 )
  8.     error(E_INTERN,"bkp_swap_elt");
  9.     if ( idx2 < 0 )
  10.     {
  11.     idx2 = r2->len;
  12.     if ( idx2 >= r2->maxlen )
  13.     {    tracecatch(sprow_xpd(r2,2*r2->maxlen+1,TYPE_SPMAT),
  14.             "bkp_swap_elt");    }
  15.  
  16.     r2->len = idx2+1;
  17.     r2->elt[idx2].col = j2;
  18.     r2->elt[idx2].val = 0.0;
  19.     /* now patch up column access path */
  20.     tmp_row = -1;    tmp_idx = j2;
  21.     chase_col(A,j2,&tmp_row,&tmp_idx,i2-1);
  22.     if ( tmp_row < 0 )
  23.     {
  24.         r2->elt[idx2].nxt_row = A->start_row[j2];
  25.         r2->elt[idx2].nxt_idx = A->start_idx[j2];
  26.         A->start_row[j2] = i2;
  27.         A->start_idx[j2] = idx2;
  28.     }
  29.     else
  30.     {
  31.         row_elt    *tmp_e;
  32.  
  33.         tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
  34.         r2->elt[idx2].nxt_row = tmp_e->nxt_row;
  35.         r2->elt[idx2].nxt_idx = tmp_e->nxt_idx;
  36.         tmp_e->nxt_row = i2;
  37.         tmp_e->nxt_idx = idx2;
  38.     }
  39.     }
  40.     else if ( r2->elt[idx2].col != j2 )
  41.     error(E_INTERN,"bkp_swap_elt");
  42.  
  43.     e1 = &(r1->elt[idx1]);    e2 = &(r2->elt[idx2]);
  44.  
  45.     tmp = e1->val;
  46.     e1->val = e2->val;
  47.     e2->val = tmp;
  48.  
  49.     return A;
  50. }
  51.  
  52. /* bkp_bump_col -- bumps row and idx to next entry in column j */
  53. row_elt    *bkp_bump_col(A, j, row, idx)
  54. SPMAT    *A;
  55. int    j, *row, *idx;
  56. {
  57.     SPROW    *r;
  58.     row_elt    *e;
  59.  
  60.     if ( *row < 0 )
  61.     {
  62.     *row = A->start_row[j];
  63.     *idx = A->start_idx[j];
  64.     }
  65.     else
  66.     {
  67.     r = &(A->row[*row]);
  68.     e = &(r->elt[*idx]);
  69.     if ( e->col != j )
  70.         error(E_INTERN,"bkp_bump_col");
  71.     *row = e->nxt_row;
  72.     *idx = e->nxt_idx;
  73.     }
  74.     if ( *row < 0 )
  75.     return (row_elt *)NULL;
  76.     else
  77.     return &(A->row[*row].elt[*idx]);
  78. }
  79.  
  80. /* bkp_interchange -- swap rows/cols i and j (symmetric pivot)
  81.     -- uses just the upper triangular part */
  82. SPMAT    *bkp_interchange(A, i1, i2)
  83. SPMAT    *A;
  84. int    i1, i2;
  85. {
  86.     int        tmp_row, tmp_idx;
  87.     int        row1, row2, idx1, idx2, tmp_row1, tmp_idx1, tmp_row2, tmp_idx2;
  88.     SPROW    *r1, *r2;
  89.     row_elt    *e1, *e2;
  90.     IVEC    *done_list = IVNULL;
  91.  
  92.     if ( ! A )
  93.     error(E_NULL,"bkp_interchange");
  94.     if ( i1 < 0 || i1 >= A->n || i2 < 0 || i2 >= A->n )
  95.     error(E_BOUNDS,"bkp_interchange");
  96.     if ( A->m != A->n )
  97.     error(E_SQUARE,"bkp_interchange");
  98.  
  99.     if ( i1 == i2 )
  100.     return A;
  101.     if ( i1 > i2 )
  102.     {    tmp_idx = i1;    i1 = i2;    i2 = tmp_idx;    }
  103.  
  104.     done_list = iv_resize(done_list,A->n);
  105.     for ( tmp_idx = 0; tmp_idx < A->n; tmp_idx++ )
  106.     done_list->ive[tmp_idx] = FALSE;
  107.     row1 = -1;        idx1 = i1;
  108.     row2 = -1;        idx2 = i2;
  109.     e1 = bkp_bump_col(A,i1,&row1,&idx1);
  110.     e2 = bkp_bump_col(A,i2,&row2,&idx2);
  111.  
  112.     while ( (row1 >= 0 && row1 < i1) || (row2 >= 0 && row2 < i1) )
  113.     /* Note: "row2 < i1" not "row2 < i2" as we must stop before the
  114.        "knee bend" */
  115.     {
  116.     if ( row1 >= 0 && row1 < i1 && ( row1 < row2 || row2 < 0 ) )
  117.     {
  118.         tmp_row1 = row1;    tmp_idx1 = idx1;
  119.         e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
  120.         if ( ! done_list->ive[row1] )
  121.         {
  122.         if ( row1 == row2 )
  123.             bkp_swap_elt(A,row1,i1,idx1,row1,i2,idx2);
  124.         else
  125.             bkp_swap_elt(A,row1,i1,idx1,row1,i2,-1);
  126.         done_list->ive[row1] = TRUE;
  127.         }
  128.         row1 = tmp_row1;    idx1 = tmp_idx1;
  129.     }
  130.     else if ( row2 >= 0 && row2 < i1 && ( row2 < row1 || row1 < 0 ) )
  131.     {
  132.         tmp_row2 = row2;    tmp_idx2 = idx2;
  133.         e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
  134.         if ( ! done_list->ive[row2] )
  135.         {
  136.         if ( row1 == row2 )
  137.             bkp_swap_elt(A,row2,i1,idx1,row2,i2,idx2);
  138.         else
  139.             bkp_swap_elt(A,row2,i1,-1,row2,i2,idx2);
  140.         done_list->ive[row2] = TRUE;
  141.         }
  142.         row2 = tmp_row2;    idx2 = tmp_idx2;
  143.     }
  144.     else if ( row1 == row2 )
  145.     {
  146.         tmp_row1 = row1;    tmp_idx1 = idx1;
  147.         e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
  148.         tmp_row2 = row2;    tmp_idx2 = idx2;
  149.         e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
  150.         if ( ! done_list->ive[row1] )
  151.         {
  152.         bkp_swap_elt(A,row1,i1,idx1,row2,i2,idx2);
  153.         done_list->ive[row1] = TRUE;
  154.         }
  155.         row1 = tmp_row1;    idx1 = tmp_idx1;
  156.         row2 = tmp_row2;    idx2 = tmp_idx2;
  157.     }
  158.     }
  159.  
  160.     /* ensure we are **past** the first knee */
  161.     while ( row2 >= 0 && row2 <= i1 )
  162.     e2 = bkp_bump_col(A,i2,&row2,&idx2);
  163.  
  164.     /* at/after 1st "knee bend" */
  165.     r1 = &(A->row[i1]);
  166.     idx1 = 0;
  167.     e1 = &(r1->elt[idx1]);
  168.     while ( row2 >= 0 && row2 < i2 )
  169.     {
  170.     /* used for update of e2 at end of loop */
  171.     tmp_row = row2;    tmp_idx = idx2;
  172.     if ( ! done_list->ive[row2] )
  173.     {
  174.         r2 = &(A->row[row2]);
  175.         bkp_bump_col(A,i2,&tmp_row,&tmp_idx);
  176.         done_list->ive[row2] = TRUE;
  177.         tmp_idx1 = unord_get_idx(r1,row2);
  178.         tracecatch(bkp_swap_elt(A,row2,i2,idx2,i1,row2,tmp_idx1),
  179.                "bkp_interchange");
  180.     }
  181.  
  182.     /* update e1 and e2 */
  183.     row2 = tmp_row;    idx2 = tmp_idx;
  184.     e2 = ( row2 >= 0 ) ? &(A->row[row2].elt[idx2]) : (row_elt *)NULL;
  185.     }
  186.  
  187.     idx1 = 0;
  188.     e1 = r1->elt;
  189.     while ( idx1 < r1->len )
  190.     {
  191.     if ( e1->col >= i2 || e1->col <= i1 )
  192.     {
  193.         idx1++;
  194.         e1++;
  195.         continue;
  196.     }
  197.     if ( ! done_list->ive[e1->col] )
  198.     {
  199.         tmp_idx2 = unord_get_idx(&(A->row[e1->col]),i2);
  200.         tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,e1->col,i2,tmp_idx2),
  201.                "bkp_interchange");
  202.         done_list->ive[e1->col] = TRUE;
  203.     }
  204.     idx1++;
  205.     e1++;
  206.     }
  207.  
  208.     /* at/after 2nd "knee bend" */
  209.     idx1 = 0;
  210.     e1 = &(r1->elt[idx1]);
  211.     r2 = &(A->row[i2]);
  212.     idx2 = 0;
  213.     e2 = &(r2->elt[idx2]);
  214.     while ( idx1 < r1->len )
  215.     {
  216.     if ( e1->col <= i2 )
  217.     {
  218.         idx1++;    e1++;
  219.         continue;
  220.     }
  221.     if ( ! done_list->ive[e1->col] )
  222.     {
  223.         tmp_idx2 = unord_get_idx(r2,e1->col);
  224.         tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,i2,e1->col,tmp_idx2),
  225.                "bkp_interchange");
  226.         done_list->ive[e1->col] = TRUE;
  227.     }
  228.     idx1++;    e1++;
  229.     }
  230.  
  231.     idx2 = 0;    e2 = r2->elt;
  232.     while ( idx2 < r2->len )
  233.     {
  234.     if ( e2->col <= i2 )
  235.     {
  236.         idx2++;    e2++;
  237.         continue;
  238.     }
  239.     if ( ! done_list->ive[e2->col] )
  240.     {
  241.         tmp_idx1 = unord_get_idx(r1,e2->col);
  242.         traceca
  243.  
  244. #ifndef ANSI_C
  245. static    void    hhldr3(x,y,z,nu1,beta,newval)
  246. double    x, y, z;
  247. Real    *nu1, *beta, *newval;
  248. #else
  249. static    void    hhldr3(double x, double y, double z,
  250.                Real *nu1, Real *beta, Real *newval)
  251. #endif
  252. {
  253.     Real    alpha;
  254.  
  255.     if ( x >= 0.0 )
  256.         alpha = sqrt(x*x+y*y+z*z);
  257.     else
  258.         alpha = -sqrt(x*x+y*y+z*z);
  259.     *nu1 = x + alpha;
  260.     *beta = 1.0/(alpha*(*nu1));
  261.     *newval = alpha;
  262. }
  263.  
  264. #ifndef ANSI_C
  265. static    void    hhldr3cols(A,k,j0,beta,nu1,nu2,nu3)
  266. MAT    *A;
  267. int    k, j0;
  268. double    beta, nu1, nu2, nu3;
  269. #else
  270. static    void    hhldr3cols(MAT *A, int k, int j0, double beta,
  271.                double nu1, double nu2, double nu3)
  272. #endif
  273. {
  274.     Real    **A_me, ip, prod;
  275.     int    j, n;
  276.  
  277.     if ( k < 0 || k+3 > A->m || j0 < 0 )
  278.         error(E_BOUNDS,"hhldr3cols");
  279.     A_me = A->me;        n = A->n;
  280.  
  281.     /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, A at 0x%lx, m = %d, n = %d\n",
  282.            __LINE__, j0, k, (long)A, A->m, A->n); */
  283.     /* printf("hhldr3cols: A (dumped) =\n");    m_dump(stdout,A); */
  284.  
  285.     for ( j = j0; j < n; j++ )
  286.     {
  287.         /*****        
  288.         ip = nu1*A_me[k][j] + nu2*A_me[k+1][j] + nu3*A_me[k+2][j];
  289.         prod = ip*beta;
  290.         A_me[k][j]   -= prod*nu1;
  291.         A_me[k+1][j] -= prod*nu2;
  292.         A_me[k+2][j] -= prod*nu3;
  293.         *****/
  294.         /* printf("hhldr3cols: j = %d\n", j); */
  295.  
  296.         ip = nu1*m_entry(A,k,j)+nu2*m_entry(A,k+1,j)+nu3*m_entry(A,k+2,j);
  297.         prod = ip*beta;
  298.         /*****
  299.         m_set_val(A,k  ,j,m_entry(A,k  ,j) - prod*nu1);
  300.         m_set_val(A,k+1,j,m_entry(A,k+1,j) - prod*nu2);
  301.         m_set_val(A,k+2,j,m_entry(A,k+2,j) - prod*nu3);
  302.         *****/
  303.         m_add_val(A,k  ,j,-prod*nu1);
  304.         m_add_val(A,k+1,j,-prod*nu2);
  305.         m_add_val(A,k+2,j,-prod*nu3);
  306.  
  307.     }
  308.     /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, m = %d, n = %d\n",
  309.            __LINE__, j0, k, A->m, A->n); */
  310.     /* putc('\n',stdout); */
  311. }
  312.  
  313. #ifndef ANSI_C
  314. static    void    hhldr3rows(A,k,i0,beta,nu1,nu2,nu3)
  315. MAT    *A;
  316. int    k, i0;
  317. double    beta, nu1, nu2, nu3;
  318. #else
  319. static    void    hhldr3rows(MAT *A, int k, int i0, double beta,
  320.                double nu1, double nu2, double nu3)
  321. #endif
  322. {
  323.     Real    **A_me, ip, prod;
  324.     int    i, m;
  325.  
  326.     /* printf("hhldr3rows:(l.%d) A at 0x%lx\n", __LINE__, (long)A); */
  327.     /* printf("hhldr3rows: k = %d\n", k); */
  328.     if ( k < 0 || k+3 > A->n )
  329.         error(E_BOUNDS,"hhldr3rows");
  330.     A_me = A->me;        m = A->m;
  331.     i0 = min(i0,m-1);
  332.  
  333.     for ( i = 0; i <= i0; i++ )
  334.     {
  335.         /****
  336.         ip = nu1*A_me[i][k] + nu2*A_me[i][k+1] + nu3*A_me[i][k+2];
  337.         prod = ip*beta;
  338.         A_me[i][k]   -= prod*nu1;
  339.         A_me[i][k+1] -= prod*nu2;
  340.         A_me[i][k+2] -= prod*nu3